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ABSTRACT 


We propose and briefly suggest how to apply the analytical tools of restricted nonlinear modal analysis 
(NMA) to problems of nuclear reactor kinetics, NPP dynamics, and NPP instrumentation and control. The 
proposed method is closely related with recent approaches by modal analysis using the reactivity matrix with 
feedback to couple neutron kinetics with thermal hydraulics in the reactor’s core. A nonlinear system of 
ordinary differential equations for mode amplitudes is obtained, projecting the dynamic equations of a model of 
NPP onto the eigenfunctions of a suitable adjoint operator. A steady state solution of the equations is taken as a 
reference, and the behaviour of transient solutions in some neighbourhood of the steady state solution is studied 
by an extension of Liapunov's First Method that enables to cope directly with the non-linear terms in the 
dynamics. In NPP dynamics these differential equations for the mode amplitudes are of polynomial type of low 
degree A few dominant modes can usually be identified. These mode amplitudes evolve almost independently 
of the other modes, more slowly and tending to slave the other mode amplitudes. Using asymptotic methods, it 
is possible to calculate a closed form analytical approximation to the response to finite amplitude perturbations 
from the given steady spatial pattern (the origin of the space of mode amplitudes). When there is finite 
amplitude instability, the method allows us to calculate the threshold amplitude as a well-defined function of 
system's parameters. This is a most significant accomplishment that the other methods cannot afford. 


1. INTRODUCTION 


Instabilities in neutron kinetics and coolant flow are important both in NPP design and NPP 
operation (Waisman, [1], Rust, [2]). In some cases, these instabilities appear after very small 
(infinitesimal) perturbations from the reactors or NPP steady state. But in other cases, the 
instabilities appear only when the perturbation moves the NPP state far enough from its 
steady state of operation. 

One of the goals of reactor design and operation is to restrict the possible states of the reactor, 


during steady operation and during transients, to remain inside a certain bounded set of 
admissible states. Also, during transients, certain restrictions must be imposed on the time 
scale of evolution of the reactor's state. So, some pertinent questions are: Will the state after a 
perturbation remain bounded within some specified set of state functions and rates of change? 
Will the state of the system return to the original steady state? Will it move farther away 
approaching to a new steady-state or to another attractor set, such as an oscillating pattern, as 
is the case of neutron-thermal-hydraulic oscillations in BWR, or xenon oscillations in large 
reactors? Or will the state of the system suffer such a severe runaway that in some instant of 
time its physical integrity will be lost? (As could happen in some unstable (at low power) 
reactor’s design, or in certain sub- critical multiplying systems). 

The best-known method to study stability problems is the First Method of Liapunov, also 
known as local linearization, or stability in the presence of infinitesimal perturbations. It 
allows us to study strictly only the stability of the linear version of the dynamic system. 
However, with the aid of a theorem due to Hartman and Grobman (Nicolis [3]), its scope can 
be extended to study the stability of non-linear systems in a small enough neighborhood of 
the steady state. But it can't cope with stability problems in the large (that is, the response to 
non-infinitesimal perturbations), like the ones posed in reactors design, operation, and 
control. 

The so called restricted analytical methods of nonlinear modal analysis (RNMA) allow us 
to tackle some of these problems. They can be considered an extension of Liapunov's First 
Method to cope directly with the non-linear terms in the dynamics. These methods 
complement both, numerical methods ab-initio, and the results obtained after applying some 
summarizing function criterion to determine a region of stability, like Liapunov's functions 
(Liapunov’s Second or Direct Method) and the like, that are often applied in theoretical 
discussions of reactors control (Stacey [4]). 


The purpose of this paper is to propose the use of certain asymptotic methods of NMA, 
mainly but not exclusively to derive thresholds of instability in nuclear power reactors. 


These analytical methods were developed to cope with finite amplitude instabilities in fluid 
mechanics. The original ideas can be found in W. Eckhaus [5] and M.Denn [6]. 

During the sixties and the seventies, the RNMA methods were applied to the study of 
stability problems in distributed parameter systems, including mass transport processes and 
chemical reactions in continuous flow chemical reactors and other systems of interest from an 
engineering standpoint. 

But RNMA methods for distributed parameter systems, such as Eckhaus’ Methods, do not 
seem to have been applied at all in the nuclear engineering field. 


It was Liapunov’s Direct Method that attracted the attention of people interested in 
mathematical methods applied to reactors kinetics. This was mainly due to the capability to 
analyze the finite —amplitude behavior of nonlinear point kinetic equations without ever 
solving them. 
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However, the main problem in practice has been the difficulties in constructing Liapunov’s 
functions which define stability regions of enough extension to be useful to introduce less 
conservative criteria for design purposes. Moreover, almost all the work was done beginning 
directly with lumped parameters models of the reactor. 


In modal analysis the different fields that give the space-time evolution of a distributed 
parameters system with a bounded space domain, are expanded in series of eigenfunctions of 
a certain suitably chosen linear operator, including the boundary conditions of the problem. 
These eigenfunctions (spatial modes) depend only of the space coordinates and are defined 
in the bounded domain. In each expansion corresponding to a given field, these spatial modes 
are weighted by unknown time dependent modal amplitudes. The modal amplitudes may be 
determined so that the series expansions represent a solution of the dynamic field equations 
verifying the initial and boundary conditions in the bounded space domain. 

The original nonlinear partial differential equations with time and space as independent 
variables are substituted by an equivalent system of nonlinear ordinary differential equations 
for the unknown mode amplitudes, with time as the only independent variable. 

From the initial conditions verified by the fields, a set of initial conditions for the mode 
amplitudes are obtained. Thus, a complex space- time field dynamics is reduced to the study 
of the evolution of a representative point in the space of mode amplitudes. 

Often it is possible to work with a relatively small number of mode amplitudes, after reducing 
the original high order dynamics to a low order one, as will be suggested below. 

To be able to apply the asymptotic methods of restricted nonlinear modal analysis 
proposed here, a proper choice of the linear operator and its eigenfunctions must be done 
starting from the nonlinear field equations for neutron dynamics. 

The critical assumption of Eckhaus’ Asymptotic Method is that the eigenvalues of this linear 
operator are widely spaced. For the nuclear reactors case the assumption is excellent, and the 
method should be widely applicable. 


Here we discuss first some background aspects of RNMA, as well as their connection with 
more common nonlinear modal methods using the reactivity matrix and the modal point 
kinetic equations. Then we give a simple example of a reactor core with a static instability of 


the “run away” type. Finally we suggest several possible applications of the proposed NMA 
methods, and we discuss some of the limitations of the methods. 


2. DESCRIPTION OF THE METHODS 


2.1. General Background 


The starting point is a nonlinear system of equations of evolution of the NPP variables, with 
focus on the reactor, and including several parameters of interest. 
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The first three equations describe the neutron kinetics in the core, using a multi-group 
diffusion approximation. Here, is the vector of neutron fluxes in each energy group. The 
variables c, are the concentration of delayed neutron emitters. And x is the state vector of 


concentration fields of relevant fission products (like xenon, iodine and samarium). The last 
two equations describe the dynamics of the thermal- hydraulics variables in the core (state 
vector w ) and the dynamics of the variables of interest outside the core (state vector w, ). 


Amongst the thermal hydraulics variables in the core, we have the temperature field in fuel, 
moderator, reflector and coolant, as well as the flow velocity fields of coolant in the core. 
Amongst the remaining variables of the NPP, we have those stemming from the primary 
circuit piping and circulating pumps, heat exchangers, turbines and electric power generating 
machines with their electric loads, or a suitable subset of these equations, depending on the 
purpose. The details of each model depend strongly of the type of NPP and of the problem 


that is going to be considered. The symbols, Bx, Ak, Xp, Xa, | and S have their usual 


meaning in multi-group diffusion theory (Stacey [4]). The scalar fission operator M and the 
vector destruction operator L are linear functions of the vector field @. Moreover, they are 


direct functions of the feedback variables x and w. Furthermore, due to the heterogeneity of 
the core, these operators are explicit functions of the position vector. 
Next we consider a steady state solution $0,¢,9,X%), Wy, W, Of the equations, after 


removing the external neutron source S. This solution corresponds to a critical state of the 
nuclear reactor core. 

The problem to be considered now is the behavior of the transient solutions in some 
neighborhood of this steady state. 


The modal methods of solution of the field equations represent the fields as linear 
combinations of known space functions weighted by unknown time functions (mode 
amplitudes). 

This tentative solution is substituted in the evolution equations. 

Applying a suitable criterion, a system of ordinary differential equations for the mode 
amplitudes is obtained. 

In the cases that we are going to consider, the set of known space functions is the numerable 
set of eigenvalues of a linear operator. 
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This operator, not necessarily self — adjoint, is constructed from the field equations. 

The criterion that is applied to obtain the evolution equations for the modal amplitudes, both 
in the common and in the restricted nonlinear modal analysis, is the projection onto each 
eigenfunction of the adjoint operator. 


If the well-known Lambda — modes are chosen as the abovementioned set of space 
functions, the equations of evolution of the corresponding mode amplitudes are a 
generalization of the point kinetics with feedback, known as modal point kinetics. 

Each mode amplitude is coupled with the others through a reactivity matrix (Turso et alter 
[7]; Ikeda et alter [8]; Ginestar et alter [9]). The feedback variables appear in the elements of 
the reactivity matrix. 


However, to apply the asymptotic RNMA methods proposed in this paper, the choice of 
the set of modal space functions must be done in a different way. 


In the equation for the neutron flux field, the following linear operator is identified, fixing 
x and w at their steady — state values: 


A,|é]= (1 8)-M[9,0, wo] [lz, — lv] £16.05 0] 


Then, the reaction diffusion equation for the neutron flux can be written as follows: 
FA "i : K 

f= Alo ]+ Ny[.x—x),0-ma]+ YA er bz +bIs 
k=l 


Here N, is anon linear operator. Fixing x — x,,w—w, , this operator is linear ing. 


It can be expanded as a sum of multi-linear operators of its arguments, beginning with a 
bilinear one. Usually a small number of multi-linear terms (two or three) will be enough. 


We pose the eigenfunction — eigenvalue problem for the operator A, : 


We use these eigenfunctions to develop the neutron fluxes o and the concentration vector of 
delayed neutron emitters c, -[v]- XZ, iM series expansions in terms of the abovementioned 


eigenfunctions. 

Then we represent the field w of feedback variables defined in the reactor core, also as a 
linear combination of suitable eigenfunctions with the corresponding time dependent mode 
amplitudes. 

Substituting this new ansatz in the evolution equations and projecting onto the 


eigenfunctions of the adjoint operator Ay, we obtain a system of nonlinear ordinary 


differential equations. 

The regularity of the non-linear operator in the original field equation has as a consequence 
that the differential equations for the mode amplitudes are of polynomial type (each second 
member is a polynomial in the mode amplitudes). 

In nuclear reactor dynamics, it is possible to work with polynomials of low degree (not 
greater than 3). 

The non-linear terms almost always involve a certain degree of coupling between mode 
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amplitudes. 

However, choosing the linear operator as suggested above, the linear term in each one 
of the equations of evolution that corresponds to mode amplitudes, appears uncoupled 
from the other mode amplitudes. The importance of this uncoupling will be seen in the 
example of reactor run away developed in 2.2. 


In order to apply an analytical approach is necessary to construct simplified mathematical 
models of the reactor core and the remaining of the NPP. It is advisable to lump parameters to 
describe the evolution of as many state variables as possible. The number of different 
“nodes” to be used depends of the frequencies of the transients that are going to be studied 
(Akcasu et alter [10]). 

In any case, we obtain a coupled system of ordinary differential equations with mode 
amplitudes as state variables. In these equations a certain number of geometric, mechanical, 
thermal, hydraulic and neutronic parameters appear. Each combination can be represented as 
a point in a parameters space. 


Depending of the problem, it is always possible to work with a finite and relatively small 
number of eigenfunctions. So, in practice, the space of mode amplitudes will be of finite 
dimension. 


To apply analytical methods, it is advisable to work in spaces of low dimensions. 

There are two general procedures to reduce a dynamic of high dimension to a dynamics of 
low dimension: the center manifold theory and the slow manifold theory. [3] 

In both cases a few dominant mode amplitudes can often be identified, so that the behavior of 
the system is essentially determined by these dominant modes. Then the mathematical 
problem may result simple enough to enable the use of the analytical tools of RNMA. 


In the center manifold theory, the stability of the equilibrium solutions is studied in the 
parameters space. The system is approximated by a linear one in each neighborhood of an 
equilibrium solution. The critical points in parameters space where the real part of one or a 
few eigenvalues change it sign are identified. In a neighborhood of these critical values of the 
parameters, central manifold theory allows us to use the corresponding modes as dominant 
modes. 

Let us suppose that the equilibrium solution is the origin of the space of mode amplitudes. 
Applying the procedure described here to choose the linear operator in the field equations of 
NPP dynamics, as in each equation of evolution the corresponding mode amplitude appears 
uncoupled from the others at least up to linear terms, the eigenvalues will be the coefficients 
of these uncoupled linear terms. If the real part of these coefficient changes its sign, the mode 
amplitude will exhibit a critical slowing down relative to the others. So, in a neighborhood of 
a critical point in parameters space, these slowed down mode amplitudes will dominate the 
dynamics. The center manifold theory can be applied to study transients near a critical state 
of the reactor core. 
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The slow manifold theory allows us to identify a few dominant modes, even far from 
criticality in parameters space. In principle it may be applied to study transients far from 
criticality in parameters space. However, to be able to apply it, a so-called slow manifold, of 
low dimension, has to be identified in the space of mode amplitudes. Furthermore, this slow 
manifold must be stable in the sense that any orbit, after a short transient, should approach 
and remain near the manifold. It is called slow because after this approach the state of the 
system changes slowly in comparison with the initial transient. The slow manifold theory 
may be applied to study transients, including nonlinear oscillations, in NPP dynamics, due to 
this happy circumstance: there is a hierarchy of widely separated time scales in nuclear 
reactor dynamics (Lewins, [11]). 

This goes from prompt neutron — dominated effects (hundredths of seconds), to heat transfer 
from fuel to coolant (tenths of seconds), coolant transit times through core and precursor 
dominated effects (seconds), coolant transit time through the entire primary circuit (tens of 
seconds), diurnal electric load variations and xenon flux tilting (tens of hours), samarium 
production (months), and fuel burn-up and transuranic isotope production (years). 

These widely separated time scales allow us to apply the methods of singular perturbation 
theory to simplify the description (Lin and Segel, [12]). 


When the variables of interest (variables of reference) belong to a certain time scale, and 
there is an attracting slow manifold, it is possible to simplify the dynamic analysis 
applying the following two principles to link scales of different orders of magnitude: 


1-The variables belonging to processes with time scales at least an order of magnitude greater 
than the reference time scale, can be considered as frozen. 


2- The variables belonging to processes that evolve with time scales at least an order of 
magnitude smaller than the variables of reference, after a short transient (produced in the so 
called inner time scale) can be considered as relaxed to equilibrium with these variables 
(evolving in the so called outer time scale). 


Using these principles, the slow manifold and a few dominant modes can usually be 
identified. These mode amplitudes evolve almost independently of the other modes, more 
slowly and tending to slave the other mode amplitudes in the same sense used in Synergetics 
(Haken, [13]). So, a high dimensional dynamic is reduced to a low dimensional one in a slow 
manifold constructed in the space of mode amplitudes. 


As we said before, to study the stability of the steady state of a nuclear reactor, we must study 
the stability of the origin in the space of dominant mode amplitudes. Using fairly well 
developed asymptotic methods (Nicolis, [3]; Denn [6]), it is possible to calculate a closed 
form analytical approximation to the response to finite amplitude perturbations from the 
given steady spatial pattern (in this case simply represented by the origin of the space of 
mode amplitudes).When there is finite amplitude instability, the method allows us to 
calculate the threshold amplitude as a well defined function of system's parameters. This is a 
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most significant accomplishment that the other methods cannot afford. 


2.2. An Example of Static Instability 


Let us consider a homogeneous bare reactor, whose extrapolated core fills a region B. This 
example is posed only to show as directly and as simply as possible how threshold 
amplitudes for instability can be derived from the methods of nonlinear modal analysis 
Because of that it is not intended to be a realistic model to be applied to some practical 
problem. For that reason, we shall work with only one group of neutrons, and with only one 
group of delayed neutron precursors. Following an approach like that employed by Tyror and 
Vaughan [14] we suppose that the neutron field flux, the precursor field and a feedback 
variable verify the equations of evolution. 


La uvigs (1 Blk + hw kw") +9, +5, (a) 
Tq 


d— $, (1b) 


T, —=Kod-w (Ic) 


Here t is time, fF is position vector, d(t,F ) is neutron flux, ¢ is a characteristic lifetime of 
neutrons, / is the delayed neutron fraction, ¢, (t,7) is a flux proportional to the density of 
delayed neutron precursor atoms, ty=//A is a characteristic time scale of precursor decay 


(with decay constant 2), M? is the quotient between the diffusion coefficient and the 
macroscopic absorption cross-section, k is the multiplication factor and S,is a source term 
due to the presence of external neutrons that are being emitted with this emission rate, w is a 
feedback variable, 7, is a characteristic time scale of evolution of w, and the parameter K 
links the neutron flux with the rate of variation of w. Furthermore, we supposed that there are 
some feedback mechanisms that make the multiplication factor k a polynomial function of w 


k=k, +kw-kw (2) 


Here ky, k; and k, are positive parameters. Consequently, if w increases from zero, k first 
increases and then decreases. The boundary conditions are d(,7,)=0 ,¢,(t,7,)=0, and 


Wt, F, ) =0 , for every ¢ and every 7, € OB (the boundary of region B). 
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Now, our purpose will be to study the stability of the stationary solutions of equations (1a), 
(1b), (1c). Let us suppose that the source of external neutrons vanishes, so 
that g, (7) =0.4,,(7)=0, and w(7)=0, is a possible steady-state solution (zero power 
solution). However, we may introduce a suitable distributed source of external neutrons, 


during a while, to produce a perturbation in the field variables relative to the strictly zero- 
power solution. After the desired perturbation is produced, the source will be removed. 


In this case the linear operator Ay [0] obtained fixing w at its steady state value w= w, is 


Ald]-“v79+(-)2¢. 


is er ( -- ewe: 


The nonlinear operator N, [¢; w] now is given by N ‘: [¢; w] = 
Then equation (la) may be cast as follows, in a form suitable to begin with NMA methods: 


~ ~ 1 1 
f= Ald]+ Ny ldsw]+—-4, +25, . 
ot Gs t 
In this case the eigenfunction- eigenvalue problem for the operator Alg| with 
homogeneous boundary conditions for ¢ is equivalent to the eigenfunction- eigenvalue 
problem for the operator -~V’ld] with the same boundary conditions forg@ . The 
eigenfunctions are the same for both operators, and the corresponding eigenvalues are related 
by a straight line. As a consequence,to study the stability of the zero-power solution, we 
represent the three fields as series of eigenfunctions of the operator —V’| | with 


homogeneous Dirichlet’s boundary conditions in the boundary of the domain B : 


b(.7)=>.4,09,@) (3a) 
w(e.7)= > B, (Ne, (7) (3b) 
#,(t.?)=YC,09.(7) (3c) 


Here—V°9,(7)= u?¢(F) , and the eigenvalues can be ordered 0 < pup <M? < WG <.... 

The eigenfunctions are orthogonal and can be _ normalized. Then 

(Qi>Pn) = | 9, (F)o,,(7)dV =6,,, being 6, Kroenecker’s 5. Substituting the ansatz (3a,b,c) 
B 
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in Eqs.(1a,b,c) and projecting onto each eigenfunction @,,, the following infinite system of 


nonlinear ordinary differential equation is obtained: 


d a n 
(A, =(- Bk 04M °), +[1- Bk, YL Ay B B, -[1- Bk, Dalek Bee, mr 


d m,n=0 m,n,q=0 d 


(4a) 


? = KA -B (4b) 


fo AP lp a thy YL dyB, hy Y Lym dyBB, |-C, (40) 


d pmng*"m~~ nq 
dt t m,n=0 m,n,q=0 


PpaON23 x. 


Here Ipmn = 1OpPmond¥ 5: 0 Tpetng =] pPmPnPqdV ands, = Jp (FSo(t,rdV. If A,(0)=0 for 


every mode index p =0,1,2,3,....and at least one of the forcing termsS, #0, then a perturbation 
will be produced from the zero-power solution. Once every projectionS, (t)has vanished, the 


perturbation is already established, and the next task is to determine its future evolution. 

Will the flux return to zero everywhere, or at least remain bounded and relatively near zero? 
Or will it suffer a severe runaway, either approaching a new steady-state solution or even 
reaching values that put in danger the physical integrity of the reactor? 

To continue, let us suppose that at zero-power the reactor is sub-critical. This means that the 


inequality k,<1+M°y, is verified, and because ww’ > uy if p #0, all the coefficients of 
the linear terms are negative in Equation (4a). Let us further suppose that 1+ M7? Lt, —k, for 
p=0 is smaller than all the other terms with p~0. Furthermore, we suppose that the outer time 
scale of the neutron flux rz, = ¢/ (1 +M? pe - ky) is an order of magnitude greater than the 


time constant of evolution of variable w and the characteristic time of precursor decay. In this 
case, applying the second principle to link time scales, we can consider that both w and gy are 
in equilibrium with the neutron flux. With this assumption, the system of equations (4a,4b,4c) 
reduces to the following: 


oA, =—(1+ Mp3 —ky)A, +R K YI AA, ~kK? YL Ay A,A, +8, (t) (5) 


pmn*~m~"n pmngq*~m*~n~"q 
t m,n=0 m,n,q=0 


Let us take as our time origin the instant in which all the external forcing terms S, vanish. 


Then, as the linear terms have negative coefficients, for a perturbation near enough to zero, 
the state of the system given by the mode amplitudes 4, (t) will return to the origin of the 


space of mode amplitudes. Now, even when we assumed that 1+M*y5 —k, is very near 


INAC 2005, Santos, SP, Brazil. 


»(0) 


zero, the coefficients of the linear parts of the equations for the other mode amplitudes are not 
near zero. This is always fulfilled in a bare core of standard dimensions. 

Let us suppose also that the external source excites mainly the zero-order mode amplitude, 
Ay (t), because S) (t) is much greater than S_, (t) for p #0. Then the amplitude A, (t) will be 


dominant, and it is possible to uncouple it from the others mode amplitudes. We obtain the 
evolution equation for the uncoupled dominant mode (This applies for t<0 if the time origin 
is taken when the perturbation is already completed) 


d 
ao = (1+. 1R ~Ky )Ay +k KT oy “kK "Togo fe S,(¢) (Zooo >0 5 Loo00 > 0) (6) 


The others mode amplitudes evolve fairly approximately this way, slaved by the dominant 
mode 


d 
— AP {1+ M72 ky JA, + KT 59942 kK Ly +S, (t) , p#0 (7) 


The others nonlinear terms are negligible (a foundation for this kind of asymptotic analysis 
can be found in references [5], [6]). 

Now, the external sources are withdrawn in the instantt=0, and we have the initial 
conditions A (0) A;(0),... with |Ay(0) much greater than the other |Ap(0) . Let us putS)(t)=0. 


The resulting homogeneous equation 


<4, = -( +M°* yu, -k,)A 9 th KT y9)A6 KK "Loo09 Ao, has at least one real root A, =0. 
t 


27 2 
Ki Loo, 000 


2~ 0000 


It is the only real root if 1+M7u —k, > . It is always stable. If it is the only real 


root, lim Ag(t)=0 and the higher mode amplitudes will be forced to approach to zero, as is 
t—+00 


easy to see from Equation (7). 


27 2 
Ki Loo 000 


2~ 0000 


But if 1+M? 5 —k, <1 there are two other real roots, Aoy andAgs. 


Ao is unstable and Agg is stable. The unstable mode amplitude is given by the formula 


K Loo = Jk Lp =A Iq (1+ MU, -k,) (8) 


Ay, = 
2KK 51 yoo 


IfA9(0)<Ao,y, then Ao(t) will return to zero dragging the other mode amplitudes to zero. 
But ifA9(0)>Ao,, Ao(t) will grow tending to Ay, and dragging the other mode amplitudes 
to follow it. We suppose that the time scales z, = ¢/ (1 + Mu - ky) of the higher order mode 


amplitudes (p40) are an order of magnitude smaller than the time scale 7, of the dominant 
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mode amplitude A(t). Then, after a short transient relative toz,, the slaved mode amplitudes 


will behave approximately as 


3S k, KT ,99-p (¢)- KK "T sooo (1) 
(4 ~ ky) 


A,(t) (9) 


So, the whole neutron flux will leave the set of states attracted by the zero flux condition and 
will tend to another stable steady-state flux. 

If the negative feedback is weak enough, the evolution equation for the dominant mode 
simplifies to 


1 A, = (1+ M72 =k, )dy +h KI qq A? (10) 
and the threshold amplitude is given by A,,, = (ent) . The stable solution Ags now has 


moved away to infinity. Given an initial condition A,(0), the solution of Equation (10) is 


given by 


(11) 


Dist) 
Herea = Ki KTow / and @ = (l4u Mo “he Y/ » SO that 2 =A,, .Then, if Ag (0)<Agy spt) 
a 


approaches zero, and for t big enough, it approaches ase. But ifAg(0)>Aoy, Ao(t) will 


runaway and for a finite time it will approach infinity. 


If we take ys as a parameter to vary, it is inversely proportional to the square of 
characteristic linear dimension of the reactors core. If the core is small enough, “) will be 
gets 
Ak T sooo 


state solution, and it will attract all the other states. So, no matter the amplitude of a 


big enough so that 1+M?u,; —k, > . Then the zero flux will be the only steady 


perturbation, it will never cause instability. 

But if the core is big enough to reverse the above inequality, a threshold amplitude appears 
and with it, an instability to finite perturbations is produced, even if the zero flux remains 
locally stable. 

The method allows us to calculate the global bifurcation of two branches, one corresponding 


to Ao, and the other to Ags, when , teaches a critical value (Figure 1). 
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Fig. 1: Bifurcation diagram for equilibrium amplitudes 


3. DISCUSSION AND CONCLUSIONS 


We gave a brief description of the procedure to obtain appropriate modal equations in order 
to apply asymptotic RNMA methods in NPP dynamics. We suggested how to reduce the 
order of the dynamics in the equations for the evolution of mode amplitudes, using the central 
manifold theory or the slow manifold theory. Then in the example of a runaway in a 
simplified model of a nuclear reactor, we saw how to find a formula for a threshold amplitude 
(when it exists) that can be used to determine the minimum amplitude of a perturbation that 
leads to instability, even when there is stability for suitably bounded perturbations. In the 
reactor model, the feedback variable w could be, for example, a difference between a local 
temperature and a reference temperature (perhaps a mean coolant temperature). In that case 
the parameter K would be given by «/(2vcF) where c_ is the specific heat capacity, F is 
coolant volumetric flow, v is the average neutron velocity in the one group description and « 
gives the heat power per neutron (see reference [11]). Then we could discuss the influence of 
coolant flow considered as a parameter, on threshold amplitudes and other aspects of modal 
dynamics in the reactor. 


With the above - now established background we can consider some problems where RNMA 
analytical methods could be of some help. This approach could be applied to derive analytical 
formulae for the stability limits and the stabilities boundaries of the fundamental and the first 
(azymuthal) mode in higher mode oscillation states of BWR. It can be applied also to Xenon 
spatial oscillations (Suarez-Antola, [15]) to derive analytical formulae for stability boundaries 
and the period of oscillation. 

The example of an idealized sub-critical reactor with a distributed external neutron source, 
developed in this paper, could be modified and extended to study the space time dynamics of 
sub-critical multiplying systems driven by external neutron sources. This could be a 
complement, from RNMA standpoint to the physics of these systems (Gandini and Salvatore, 


[16]). 


The methods like Liapunov's Direct Method, allow the construction of a priori bounds for 
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stability without a detailed consideration of the evolution of the state of the system. On the 
contrary, RNMA focuses as much as possible in a description of the dynamics. Instead of 
identifying a region around the steady state such that any state located there will be attracted 
to the steady state, or at least will remain in the region, RNMA gives estimates of the onset of 
instability: it allows the identification of the perturbation amplitudes that take out the state of 
the system from the admissible region. This could have practical importance in relation with 
reactor design and operation, because one of the design and operation goals, as already said, 
is to maintain the state of the reactor inside a set of admissible states. Now, regions 
determined by other methods may be unduly restrictive, in comparison with the regions 
determined applying RNMA. 

Due to pragmatic and understandable reasons, the design criteria usually adopted in many 
countries are, even know, conservative. Most of them guarantee stability under almost all 
conceivable operation conditions. So that, in at least the first stages of design often linear 
stability tools are enough, perhaps except for early safety studies. 

The price that this has is the higher costs of construction and operation, so less conservative 
criteria are being developed. In this process, RNMA could be of some help. 

Last, but not least, we can mention some of the main limitations of NMA methods. As 
Morton Denn [6] says, these are asymptotic methods that depend of some sort of regularity in 
the non-linear operator of the field equation. This regularity produces some kind of continuity 
of behavior between the linear and the non-linear regimes. The danger of this is that if the 
phenomena of interest are produced beyond the range of validity of the approximations, they 
could pass unnoticed. 

The other disadvantage is that if there aren't a few dominant modes the analytical approach 
probably will not succeed, and a numerical calculation should be undertaken. 

In this point it is necessary to assess if it is better to leave RNMA equations and use a 
computer code for an ab-initio discretization of field equations. 

In any case, RNMA methods that allow an estimation of the response to finite-amplitude 
disturbances using an extension of Liapunov’s First Method are straightforward in principle. 
They can provide a considerable amount of information, of value by itself and for the design 
of digital simulations applying computer codes corresponding to more complex and realistic 
mathematical models of the nuclear reactors. 


The calculations needed to apply the asymptotic methods of RNMA may be fairly long and 
tedious. However, powerful symbol manipulation packages are now available to develop 
complex symbolic calculations in the computer. 

Of course, these tools were unavailable when the original approach to RNMA was 
constructed. Therefore, daunting pencil and paper calculations had to be undertaken during 
the sixties and the seventies to apply RNMA methods. Then, it was natural that people shifted 
to other analytical tools easier to use, and to numerical methods, applied from the very 
beginning. 

Nevertheless, new and useful results can be obtained making digital simulations with the 
ordinary differential equations obtained from RNMA of partial differential equations. A 
famous example is the numerical discovery of chaos, at the beginning of the sixties, during 
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numerical simulations of the dynamics corresponding to Lorentz equations. These nonlinear 
ordinary differential equations are the equations of evolution of the coupled mode amplitudes 
obtained from a very simplified model of an atmosphere. 


The situation now is different. Even if present day symbolic packages may be still crude for 
certain specific manipulations, the processes of improving them are going on very fast. 
Furthermore, the contemporary tools of nonlinear science and a very significant amount of 
experience in applying them is already available. 

So, it seems that the proposal to apply RNMA methods to study stability problems in NPP, 
made in this paper, may be timely enough. 


— 
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Synopsis: It was proposed recently to apply an extension of Liapunov’s first method 
to the non-linear regime, known as non-linear modal analysis (NMA), to the study 
of space-time problems in nuclear reactor kinetics, nuclear power plant dynamics and 
nuclear power plant instrumentation and control. The present communication suggests 
how to apply NMA methods to the study of Xenon spatial oscillations in large nuclear 
reactors. The set of non-linear modal equations derived by J. Lewins for neutron flux, 
Xenon concentration and Iodine concentration are discussed, and a modified version of 
these equations is taken as a starting point. Using the methods of singular perturbation 
theory a slow manifold is constructed in the space of mode amplitudes. This allows the 
reduction of the original high dimensional dynamics to a low dimensional one. The 
equation of evolution that correspond to the amplitudes in the fundamental mode are 
decoupled from the first harmonic. An analytical formula is obtained for the case of 
gross (in phase) xenon oscillations. It gives an instability threshold for finite 
perturbations from a stable steady state. Then, the fundamental mode amplitudes are 
fixed at their steady state values. An analytical formula is obtained from the equations 
of evolution of the amplitudes corresponding to the first harmonic. It gives an instability 
threshold for out-of-phase (spatial) xenon oscillations. Some suggestions are given in 
order to derive better and deeper analytical and numerical results from the whole set of 
non-linear evolution equation posed in this paper. These results could be applied to the 
discussion of neutron flux and temperature excursions in critical locations of the 
reactor’s core. The results of NMA can be validated from the results obtained applying 
suitable computer codes, using homogenization theory to link the complex 
heterogeneous model of the codes with the simplified mathematical model used for 
NMA. 


Key words: xenon gross and spatial oscillations, nonlinear modal analysis, nuclear 
reactors dynamics 


1. Introduction 


Instability problems in nuclear power reactors pose important restrictions in design and 
operation of nuclear power plants (NPP). Recently it was proposed to apply an 
extension of Liapunov’s First Method, known as nonlinear modal analysis (NMA), 
to study stability problems of nuclear reactor kinetics, NPP dynamics, and NPP 


instrumentation and control [1]. The purpose of the present communication is to suggest 
how to apply NMA to the study of xenon oscillations. 


The effects of Xenon are a potential source of instability in thermal reactors. 

Xenon oscillations are produced by the delay between xenon burn-up and xenon build 
up from iodine decay when a change in local neutron flux produces an imbalance 
between both processes. This can establish an oscillatory regime in reactor power with 
dangerous peak values, mainly when the reactor is used for load following, with 
frequent power changes. 


There are two types of oscillations: the so called in-phase or gross xenon oscillations, 
and the so called out-of- phase oscillations. 

During gross oscillations the reactor power changes in phase in every point of the core. 
In out-of-phase oscillations, local xenon concentrations and power go up in one region 
of the core and down in another region. 

Control systems are designed to limit the amplitude of both types of oscillation and to 
stabilize the state of the reactor near a chosen steady state. 


Since the discovery of the problems posed by xenon instabilities in nuclear power 
reactors, a lot of work has been done about them, using both analytical and numerical 
methods. 

Gross xenon oscillations were the first to be studied using point reactor kinetics with 
temperature and xenon feedback. Most of the early work applied linear stability analysis 
coupled with digital simulations of the dynamics [2]. 

During the sixties and the seventies, some work was done applying the asymptotic 
methods of averaging to gross xenon oscillations [3]. This method, originally developed 
in the former Soviet Union, allows an analytical approach to the nonlinearities in the 
point kinetics equations with feedback, in the same spirit as the NMA methods 
proposed in the present communication and in reference [1]. However, the calculations 
were very complicated and somewhat clumsy at that time, so that this analytical 
approach to gross xenon instabilities was interrupted. 


As gross xenon oscillations are slow, it was relatively easy to control them in steady 
state operation. However, the need to follow the electric load as closely as possible, 
poses xenon related stability problems during power transients. 

In case of a complete loss of the load, if the NPP has the capability to bypass a 
significant part of its full-power steam load, it is advisable to be able to execute a partial 
trip in order to reduce the time to return to full power operation. 

The control of the reactor in these transient conditions requires a real-time knowledge of 
xenon concentrations. As consequence, research using nonlinear point kinetics models 
and sophisticated computer algorithms for system identifications is been done even 
today [4]. 


The need to design and to construct larger nuclear power reactors brought the problem 
of the out-of-phase or spatial xenon oscillations. The first analytical approach to spatial 
instabilities in NPP was done using linear modal analysis [5]. The stability of the first 
uncoupled spatial modes was studied using the methods of linear control theory and 
applied to several types of thermal reactors: graphite moderated and gas cooled, light 
and heavy water moderated and light and heavy water cooled. 

In a cylindrical core, xenon spatial oscillations may be in principle of three types: radial, 
azimuthal and axial. While radial and azimuthal oscillatory modes are not so easily 
excited, the axial oscillations are often easier to excite. As consequence of the 


importance of these oscillations in PWR, most of the work was done in relation with 
axial oscillations [6] [7] [8]. 

Onega and Kisner developed a two-point axial xenon oscillation model using one group 
diffusion approach and made digital simulations of the dynamics [6]. 

Song and Cho constructed a two-group and two-point axial xenon oscillation model [7]. 
These two-point models split the reactors core in two regions and the neutron flux field 
and the xenon and iodine concentration fields are averaged over each region. 

After averaging, a system of coupled nonlinear ordinary differential equations is 
obtained to describe the evolution of the amplitudes of the fundamental and the first 
spatial modes. Jointly with computer codes, these models are being used to design 
control systems in PWR [8]. 

The NMA approach to spatial xenon oscillations proposed in this paper, is closely 
related to modal kinetic analysis using static Lambda-modes [9] and modal reactivities 
[10] [11]. It does not imply spatial averaging over two reactor regions as the above 
mentioned two-point kinetics models. 


2. A framework for nonlinear modal analysis of xenon instabilities 
2.1. A summary of NMA methodology to study xenon effects 


To suggest how to apply NMA methods to xenon instability problems, the one energy 
group approach to neutron kinetics may be enough. We use the following cae 
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Time, as usual, is represented by ¢ and the position vector by /. In the first equation, is 
the scalar neutron flux, c, are the concentrations of delayed neutron precursors, and £ 
is the fraction of delayed fission neutrons. In the equations of evolution of the 
precursors, /, are the group fractions of delayed fission neutrons, and A, are the decay 
constants of precursors. In the balance equations for iodine and xenon, x is the 
concentration of xenon, and i is the concentration of iodine. g, and @, are the iodine 
and direct xenon yields and 4, ,A, are the corresponding decay constants. In the 
equation that gives the evolution of the local temperature AJ, 7 is temperature’s time 
constant, that represents the time lag in the temperature feedback produced by the 
thermal capacity of the reactor. A7, and ¢, are a reference temperature and neutron 


flux respectively [5]. M [o| =v) ,¢@ is the production operator, being v the average 
number of neutrons per fission and », the macroscopic fission cross-section. We 


neglect the feedback effects in this last cross-section. 
LI¢, XL. | = om [o|- o ,x@—QAT®@ is the destruction operator, where o, is xenon’s 


microscopic absorption cross-section and a@ is a temperature feedback coefficient. 
i [d| =-Ve (DV ¢)+ ~, ¢, being D the diffusion coefficient and >’, the macroscopic 


absorption cross-section in absence of xenon and in zero power conditions. As we are 
going to describe flux variations in a time scale of hours, the description of the reactor 
can be simplified if we treat the delayed neutrons in a quasi-static approximation. Thus 
it is possible to derive, using a procedure explained in reference [9] (page 215): 
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Here rz, is the average life time of neutron precursors. To simplify the calculations we 
consider now a homogeneous core so that all the parameters are space independent. 
Besides we consider the case of **U where 8 =0.0065, 7, =12.85, so that 
A, = ft, =8.3x10~s is an effective generation time, and the prompt neutron 

1 


UV f 


generation time is A, = =10~“s. In this case the neutron evolution equation 


can be further simplified to A.M |] = M,|¢|-L¢,x,AT]. This last equation may be re- 


at 
written thus: 


ap _ Mo 


V°d+plx, AT 1 
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Here M) == is the neutron migration area, k, = = is the zero-power 


multiplication constant for an infinite medium. 


An effective reactivity with feedback (being p,, =1-— x? is given by the following 


Ox AT 
formula: plx, AT |= ?.-—=- = (2) 
voy Voy 
2 
In —* we can substitute k, by 1, because we are considering a reactor that is 


subcritical in relation with prompt neutrons. Then to study xenon instability problems 
we can start from the following equations: 


A, = =M>V°¢+ plx, AT (3a) 
Oi 
Ot =@; 2g p-A,i (3b) 
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The reactivity plx, AT | is given by equation (2). If B is the region of the extrapolated 
core of the reactor, ¢ must be zero when r belongs to the boundary of B. 


The problem now is to study the stability of the steady-state solutions of these dynamic 
equations. Given the homogeneous boundary condition applied to the neutron flux, the 
zero solution is always an equilibrium solution. But if p, is positive and the 


dimensions of the core are big enough, there is always a unique positive solution ¢, (r), 
x: (r), i (r) and AT, (r) that represents the corresponding just-critical state of the 


reactor. The neutron evolution equation is posed as follows: 


A. f= Ayld}+ No[tsx-,AT-A7,] (4) 


A,|¢]=M2V°¢+ plx,,AT, pis a linear operator (self-adjoint). 


N,[¢;x—x),AT - AT, ]= " (x-x,)p = (AT-AT,)é is a nonlinear one. 
f f 


The idea now is to consider the eigenfunctions ‘¥,, and the eigenvalues @, of the three- 


dimensional Sturm-Liouville problem Ay [w]= o-Y with Y defined in B and with 


zero boundary conditions. 
The numerable set of eigenfunctions is rich enough to allow us to represent the neutron 
flux, the concentrations of xenon and iodine, and the local temperature, as follows: 
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Besides the eigenfunctions are orthogonal in the sense that | w (rye r)d )dV = 0 if 
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m#n and can be normalized. 
Substituting the ansatz (5) in the dynamic equations (3) and projecting onto each 
eigenfunction Y , (r), an infinite set of nonlinear ordinary differential equations for the 


time dependent mode amplitudes is obtained. These equations are a modification and a 
remarkable generalization of the equations derived by Lewins [12] (pp 111-114) in his 
approach to spatial xenon oscillations. 


The operator Ay [¢] was chosen such that the mode amplitudes ® i (t) for the neutron 


flux appear uncoupled to the first order (that is, uncoupled to linear terms). The reasons 
for doing this are explained in [1]. 

In the modal expansion, the eigenfunctions are ordered according to the magnitude of 
their respective eigenvalues. Then, the expansion is truncated at a certain 
eigenfunction Y,, (r), using the following criterion. The higher harmonics will be 


excited only by reactivity perturbations localized in a region small enough, so that the 
reciprocal of a representative dimension of this region is longer than the eigenvalue 
separation of these harmonics. As consequence, modal expansions with only a few 
eigenfunctions may be almost as accurate as the output of full three-dimensional modal 
methods, if the initial reactivity perturbation is not too localized. For this reason, low 
order modal expansions may be enough to study common xenon instability problems 
but are no adequate to describe a transient produced by, for example, a sudden drop of a 
control bar. After truncation above the order M, we obtain a finite system of nonlinear 
ordinary differential equations to describe the evolution of the mode amplitudes in a 
finite dimensional euclidean space. To study the dynamics in this state space, besides 
numerical studies, we have at our disposal the powerful set of analytical tools of 
nonlinear science, including asymptotic methods of analysis [13]. A combination of 
these last methods with the center manifold theorem or with the so-called slow manifold 
theorems, allows a further reduction to lower dimensional dynamics [14] and the 
derivation of formulae for thresholds of instability. Some examples are given in [1] and 
here in section 2.3. 


2.2. Nonlinear modal equations for xenon oscillations 

From now on, we take the zero solution (zero power state) as reference. Of course, this 
is not the best thing to do in order to study instabilities during the operation of the 
reactor. However, it has two advantages. First, we obtain a set of equations that, in the 
case of gross xenon oscillations can be easily compared with the equations already 
studied by others. Second, because we can see how it is possible to obtain another 


steady state solution, different from the origin, and study their stability using the same 
modal framework for all of them. To fix ideas, let us consider the case of a PWR, 
withA, ~0.ls, and 7 3s. This time scale of temperature is much smaller than the 


time scale of iodine (7hs), and the local time scale of xenon. The local time scale of the 
neutron flux, in spite of the value of the effective generation time, may be much larger 
than . 


M,V’ 
It is of the order of the time a/[ - Z + plx, ar) 


If the state of the reactor is near an attainable steady-state, the denominator may be 
small because both the buckling and the effective reactivity have opposite signs and 
almost the same absolute values. However, if the state of the reactor is far enough from 
a steady-state, the time scale of evolution of the flux may be small. 

In any case, we can consider that the local temperature is relaxed to its equilibrium 


AT, 
value AT, =——-¢. 
aa’, ; Sant 
Introducing the new parameter y= ro and applying the procedure described in 
VD Px 


section 2.1, and taking the first two eigenfunctions only, we obtain for the mode 
amplitudes: 


For the fundamental mode: 
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Notice the coupling between the fundamental and the first harmonic in the equations of 
evolution for the modal amplitudes corresponding to neutron flux and xenon 
concentration. 


For the first harmonic: 
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Observe the nonlinear coupling between the fundamental and the first harmonic in the 
equations of evolution for the mode amplitudes corresponding to neutron flux and 


xenon concentration. In this equations p, = p,, -M> uo and p* =p, —M>u,; , where 
HM, and yu; are the fundamental and the first harmonic eigenvalues of —V*Y = wu? 
with yw =0on the boundary. 


The relation between one of this eigenvalues, ° , and the corresponding eigenvalue @, 
of the operator Ay introduced in section 2.1 above, is @, =p? =p,,-Mou. 
p, is the so called “static reactivity” of mode n, without feedback. By definition 
Pm =[¥,(r)¥,, (7) ¥,, (1) aV for p,m,n =0,1 
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From its definition we see that @,,,,, remains invariant under any permutation of its sub- 


indexes. 
Y, is symmetric and positive and , is anti-symmetric considering a mid plane 


through the core, so that Ao) >9, Ao, =O, Ai; > 9, A), =90. 


2.3. Thresholds of instability for oscillations 

Using only two eigenfunctions to approximate the fields in the reactor’s core, we 
represent the neutron flux by ®,(t)¥,(r)+@,(¢)¥,(r) and iodine’s and xenon’s 
concentrations by i,(t)¥,(r)+i,(¢),(r) and x,(¢)¥,(7)+~x,(t)¥;,(r) respectively. At 
this level of description, in-phase xenon oscillations appear mainly in relation with the 


amplitudes of the fundamental mode, ®,(¢),x,(¢),i)(¢) , while out-of-phase oscillations 
appear mainly related with the amplitudes of the first harmonic ®, (¢), x, (0), A (t). As a 
first approximation we may assume that in gross oscillations only the first mode is 


excited, and in out-of-phase oscillations the fundamental mode remains near its steady 
state and only the first harmonic is excited. 


2.3.1 Gross xenon oscillations: an instability threshold to finite amplitude 
perturbations 

If we neglect the nonlinear coupling terms in equations (7a) and (7c), we obtain a 
simplified set of equations. This set is identical, from a mathematical point of view, 
with the set of equations used in early analytical and numerical studies of gross xenon 
oscillations [2] [3]. As consequence all the results already obtained can be translated to 
the dynamics of the uncoupled fundamental mode. If the state of the system is not too 
far from a steady state, it is possible to suppose that ®, (t) is always relaxed to 


equilibrium with x, (t). Then, applying a procedure analogous as the one used in [2] it is 
possible to derive the following equation of evolution: 
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Here by definition x, (t)= x,(1+€(t)), where ¥, is the non-trivial steady state value of 


the modal amplitude of xenon concentration. The parameters in equation (9) are defined 
as follows: 
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®, is the non-trivial steady state value of the modal amplitude of the neutron flux. It is 
an increasing function of the modal reactivity 5. 

Both @; and bc are always positive. However, the linear damping coefficient a may 
change its sign. 


Let us suppose that it is positive, so that the steady state of the reactor is stable to 
infinitesimal perturbations (in this case the state of the reactor is given simply by the 
variables &(t) ang 2) ). 
t 

Let us further suppose that the parameter a is small enough to enable the application of 
asymptotic methods. If time ¢ is substituted by —t in (9), the equation that results is 
mathematically identical to an equation posed by Bautin and studied by multiple time 
scales methods in [13] (pp. 402-406). After an inversion in the direction of time, an 
unstable rest point becomes stable, and a stable limit cycle becomes unstable. So, from 
the results derived in [13] we obtain for equation (9): in the plane of the variables E(t) 


dé(t) 


and ae. the origin is stable and there is an unstable limit cycle with center in the 
t 


origin and radius 2Va/bce . 

The origin is the steady state of the reactor at this level of description. 

If a perturbation from the steady state leaves the state of the reactor inside the circle, it 
will return to the origin. But if the perturbation leaves the state of the reactor in the 
region outside the circle, the state will move farther away. So, a subcritical Hopf 
bifurcation type of instability is produced. 


2.3.2 Out-of-phase xenon oscillations: instability threshold and period 
Let us consider now the excitation of the first harmonic while the modal amplitudes of 
the fundamental mode remain at their steady-state values. Then, equations (8) become: 
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As we did in the case of gross xenon oscillations, we will assume that the modal 
amplitude of the neutron flux is relaxed to equilibrium with the modal amplitude of 
xenon concentration: 


0,04 [,S aus} /[A9( $8 +. |-o x,(t) (11) 
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We thus obtain a system of linear differential equations in the mode amplitudes 
corresponding to xenon and iodine concentrations. Following a procedure similar to that 
employed to derive equation (9), we find the equation for a damped harmonic oscillator 
2 
zd Plea =0 
dt dt 
The damping coefficient a is given by 
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This coefficient may change its sign if the temperature feedback is weak enough and the 
flux is high enough. 

The combinations of values of the parameters of the reactor such that a is zero defines a 
stability boundary for infinitesimal perturbations. The natural frequency is given by: 
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The true frequency of the oscillator and its period may be calculated by well known 
formulae for a damped linear oscillator. 


3. Discussion and conclusions 


In this communication we used only two eigenfunctions to represent the spatial 
distribution of the fields in the reactors core. Without any assumption related with the 
slow manifolds, the dimension of the state space corresponding to the mode amplitudes 
is eight. When the local temperature is relaxed to equilibrium with neutron flux, this is 
reflected by two linear relations, one between the amplitudes for the fundamental mode, 
and the other between the corresponding ones for the first harmonic mode. Thus, the 
dimension of the manifold in which the dynamic behavior of the system is produced, is 
reduced to six. This makes sense because, after a short transient, the state of the system 
approaches to and remains near the lower dimensional manifold, where it evolves much 
more slowly. Not every slow manifold is stable in this sense, so that this issue must be 
studied carefully. 

In the first part of section 2.3, we uncoupled the fundamental mode from the first 
harmonic neglecting the coupling terms. Thus, we reduced the dimension of the state 
space from six to three. This allowed us to consider in-phase xenon oscillations in a 
mathematical framework analogous to the framework obtained applying the equations 
of point kinetics with feedbacks. However, in our case, we have a well-defined 
connection between the static reactivity and the geometry of the core. Besides, we 
obtained an analytical formula for a threshold of instability to finite perturbations that 
relates the amplitude of a threshold perturbation with the parameters of the reactor. The 
consequence of this could be worth to be studied. 

In the second part of section 2.3, we uncoupled, in an indirect way, the first harmonic 
from the fundamental mode. Instead of neglecting the coupling terms, we substituted 
the amplitudes corresponding to the fundamental mode by their steady-state values. A 
three-dimensional state space resulted from this. A further reduction in the dimension of 
the state space was obtained by the assumption that the flux amplitude and xenon 
concentration amplitude are in equilibrium. This gave us a two-dimensional linear 
system equivalent to a damped harmonic oscillator. The damping coefficient and natural 
frequency was given by a well-defined analytical relationship in terms of the parameters 
of the reactor. These results are a generalization in the framework of modal analysis, of 
the results for xenon spatial oscillations obtained using two-nodes nodal models [15]. 
The linear approximation to the dynamics of the first harmonic mode is not satisfactory. 
The full nonlinear system should be studied, without neglecting the influence of small 
oscillations in the amplitudes corresponding to the fundamental mode. An effect related 
with the well known “parametric resonance” could be produced. 
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ABSTRACT 


The cores of nuclear reactors, including its structural parts and cooling fluids, are complex mechanical systems able to 
vibrate in a set of normal modes and frequencies, if suitable perturbed. The cyclic variations in the strain state of the 
core materials may modify the reactivity, and thus thermal power, producing variations in strain due to thermal-elastic 
effects. If the variation of the temperature field is fast enough and if the Doppler Effect and other stabilizing prompt 
effects in the fuel are weak enough, a fast oscillatory instability could be produced, coupled with mechanical vibrations 
of small enough amplitude that they will not be excluded by the procedures of conventional mechanical design. After a 
careful discussion of the time scales of neutron kinetics, thermal-elastic and vibration phenomena, a simple lumped 
parameter mathematical model is constructed in order to study, in a first approximation, the stability of the reactor. An 
integro-differential equation for power kinetics is derived. Under certain conditions, fast oscillatory instabilities are 
found when power is greater than a threshold value, and the delay in the global power feedback loop is big enough. 
Approximate analytical formulae are given for the power threshold, critical delay and power oscillation frequency. It is 
shown that if prompt stabilizing fuel effects are strong enough, dangerous fast power oscillations due to mechanical- 
thermal-nuclear coupling phenomena can not appear at any power level. 


1. INTRODUCTION 


The cores of nuclear reactors, considered as mechanical systems, are able to vibrate in a set of 
normal modes and frequencies, if suitable perturbed. The time variations in the strain state of the 
core materials may produce changes in density. Changes in density modify the reactivity. Changes 
in reactivity modify thermal power. Modifications in thermal power produce variations in 
temperature fields. Variations in temperature produce variations in strain due to thermal-elastic 
effects. If these strain variations are fast enough, inertial effects must be taken into account, and 
mechanical vibrations of the core’s materials may be produced. 

Thus the cycle of events may be closed. However, under normal conditions, in a suitably designed 
core, the global effect of this feedback tends to increase the stability of the reactor. 

But if the variation of the temperature field is fast enough and if the Doppler Effect and other 
stabilizing prompt effects in the fuel are weak enough, a fast oscillatory instability could be 


produced, coupled with mechanical vibrations of small enough amplitude that they will not be 
excluded by the procedures of conventional mechanical design. 

It seems that this kind of oscillatory instability, with frequencies in the scale of tenths of seconds or 
less, could have been part of the processes that occurred during several extreme experiments (like 
the series named with the acronyms BORAX and SPERT and done with research reactors), during 
certain incidents (like certain power oscillations sometimes found in high temperature gas cooled 
reactors of old design) or even during true accidents, like certain reactivity accidents in both 
research and power reactors. 

In reference [1] a simple lumped parameter mathematical model (a nonlinear point kinetics model 
with several feedbacks) was constructed in order to study the stability of a reactor. 

It may be considered as a generalization of a model due to Thompson that appeared in 1988 [2]. 
Thompson’s model is at its turn a generalization, to introduce thermal-mechanical coupling, of a 
classical model proposed by Weinberg and Wigner in reference [3] to describe short (compared 
with the period of the delayed neutrons) power excursions. The generalization introduced in 
reference [1] takes into account certain prompt feedback effects not considered in previous models 
related in a way or another with some kind of thermal-mechanical coupling, like the 
abovementioned Thompson’s model. 


In the present paper, an equivalent single nonlinear integro-differential equation for power kinetics 
is derived, with a continuous distribution of time lags. 

Then this equation is approximated by a nonlinear equation with a discrete lag. 

This last equation may be considered as a major simplification of the two coupled oscillator model 
introduced in reference [1]. Nevertheless, it retains some of its main predictive capabilities. 

Under certain conditions, fast oscillatory instabilities are found when power is greater than a 
threshold value, and the delay in the global power feedback loop is big enough. 

Approximate analytical formulae are given for the abovementioned power threshold, critical delay 
and power oscillation frequency. 

In the framework of the mathematical model, it is shown that if prompt stabilizing fuel effects are 
strong enough, dangerous fast power oscillations due to mechanical-thermal-nuclear coupling 
phenomena can not appear at any power level. 


2. THERMAL-MECHANICAL COUPLING AND FAST OSCILLATORY 
INSTABILITIES IN NUCLEAR REACTORS 


There is a hierarchy of widely separated time scales in nuclear reactor dynamics. This goes 
from very fast transients related with local effects of nuclear fission in fuel materials (pehaps 
microseconds, [4]), to prompt neutron — dominated effects (hundredths of seconds), heat 
transfer from fuel to coolant (tenths of seconds), coolant transit times through core and 
precursor dominated effects (seconds), coolant transit time through the entire primary circuit 
(tens of seconds), diurnal electric load variations and xenon flux tilting (tens of hours), 
samarium production (months), and fuel burn-up and transuranic isotope production (years). 


When the variables of interest (variables of reference) belong to a certain time scale, it is 
possible to simplify the dynamic analysis applying the following two principles to link scales 
of different orders of magnitude: 

1-The variables belonging to processes with time scales at least an order of magnitude 
greater than the reference time scale, can be considered as frozen. 

2- The variables belonging to processes that evolve with time scales at least an order of 
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magnitude smaller than the variables of reference, after a short transient (produced in the so 
called inner time scale in the terminology of singular perturbation theory [5]) can be 
considered as relaxed to equilibrium with these variables (evolving in the so called outer time 
scale). 


Let us begin with a reactor, generating a steady thermal power P, and with the remaining state 
variables at their steady sate values. Then at the time origin a sudden perturbation is produced 
in the initial power. The time scale of reference, for the transients we are interested in is of 
the order of hundredths of seconds in a thermal reactor, so in the equation for power point- 


kinetics with feedback the term due to delayed neutrons appear as a steady source A. P, 


being / the fraction of delayed neutrons and A the mean time between neutron generations. 
The reactivity o will be a function of the temperatures and densities through the 


corresponding feedback coefficients. However, some temperature effects on reactivity are 
very fast in comparision with our time scale of reference (this last scale is of the same 


: Dee 4s : : : : 
numerical order as —), like certain prompt effects in uranium oxide fuels due to an almost 


instantaneous Doppler broadening in the central region of the relatively poor thermal 
conducting pellets of uranium oxide. The lag between power and the temperatures related 
with prompt effects may be neglected, so that these temperatures can be considered as 
relaxed to equilibrium with the instantaneous power. Temperatures like average fuel 
temperature evolve near the reference time scale, and the lag between them and the 
instantaneous power can not be neglected. Others, like average coolant temperatures, evolve 
in scales of much higher order, so that these temperatures may be considered as frozen and 
the heat removal by the coolant may be taken as a constant rate P, during the short power 
excursions studied here. 

The coupling between variations in densities and variations in temperatures necessarily must 
include inertial effects if the time scales of mechanical vibrations in the core structures, t,, , 1S 


at least of the same order of magnitude as both the time scalest, of power variations and the 


characteristic thermal time scalest,of those same structures[6]. If G ~ 0.007 and A = 10s, 
: me we te 2 
then the time scale of thermal power variations is — ~ 10° seconds for a thermal neutrons 


reactor. (It is a thousand times smaller for a fast neutrons reactor). Flexural vibration modes 
of structural solids like slender bars and thin plates, as well as certain acoustic modes in two 
phase fluids like water with vapor bubbles, may have frequencies low enough to comply with 
the requirement mentioned above in thermal neutron reactors. 


2.1. Generalization of the Weinberg-Wigner-Thompson Model of Point Kinetics 

Let us describe the reactor using the power P(t), a representative average temperature 7 (rt) 
and a representative average density d(t) of the core. Then, taking into account the precedent 
discussion about time scales we pose, being P, ,7, and d, the steady state values: 


ages ae eee 
ate p)Ree iz (1) 
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d(t)—d 
If y(t) = au , we suppose that the reactivity is given by (6p, is the reactivity added 


0 
from outside, by control bars movement or any other external mechanism): 


p=p,t+a,(P—P,)+a,(T-T,)+a,-y (2) 


The feedback reactivity coefficients are: a@,(prompt power) a, (thermal) a, (density). In 
order to have static stability we assume that both a, anda, are negative and that @ , is positive. 


For a constant heat removal rate P,, ifC is a suitable thermal capacity, the representative 


temperature is given by: 
Cc: ae P-P (3) 
dt : 


The coupling between density and temperature, including inertial effects is given by [1], [2]: 


2 
- boy 2 +0, y == 0," BAT ~Ty) (4) 
Here b is a thermal expansion coefficient, c,, =2-¢,,, °@,, 18 a mechanical damping parameter 
(being ¢,, the mechanical damping ratio) and @,, is a natural frequency of mechanical 
oscillations. Ifa@,=0,a,=Oanda, #0we obtain Weinberg and Wigner’s model. 
Ifa, #Oanda, #0we have Thompson’s model. Witha, #0,a, #0 anda, # Owe obtain 
the generalization proposed in this paper. All the parameters will be considered as constant. 


Now let us introduce the new variables v = < and x = mf 2 and define: 
t 0 


The new functions: f(x)=e,.e7% +c,.e** (5a) g(x)=e" -1 (5b) 

The nuclear damping parameters: Cy = f (6a) Cr= a (6b) 
= 0700 

The thermal frequency: Or mye a) (6c) 


A.C 
Eliminating the temperature from equations (1) to (4), and taking into account definitions 
(5a), (Sb), (6a), (6b) and (6c), we derive the following second order differential equations: 


d° d. a, 
5 rt f(x) 5 -+0,".¢(x)= rn v (7a) 
d°y dv o,, b.P, 

eae =m 7b 
dt’? =” dt TOn C (x) m) 


A nonlinear nuclear-thermal oscillator given by (7a) (that corresponds to reactor point 
kinetics with thermal-elastic feedback and with frozen delayed neutron effects) is coupled 
nonlinearly with a linear mechanical-thermal oscillator given by (7b) (that corresponds to 
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the first normal mode of suitable mechanical vibrations excited by thermo-elastic effects). 


The steady state of the reactor corresponds now to x =0, < =Oandv=0, < =0. 
t t 


2.2. A Nonlinear Integro-Differential Equation for Fast Power Transients 
Taking the Laplace Transform in both members of equation (7b), with zero initial conditions 


for v and solving the resultant equation for the Laplace Transform v(s)of v(t), and 
t 


inverting the transform V(s), we obtain: 


* DP, | 
v(t) = aoe 0 | G(t —t')g(x('))- de (8) 
1 At Ayt 
G(t)= GA) .(e e ) (9) 
Ing = Oy (=n t VE? -1) (10) 


(If¢,, =1, thenG(t)=t-e"). 
Substituting (8) in (7a), we derive the following nonlinear integrodifferential equation, 


being @, a density related oscillation frequency and K (t)defined such that } K (t) -dt=1: 
0 


re FS + 0,?.0(0)4 Oy J EGCEECU) ars0 (11) 
oO, Geely (12) 
: C-A 


Whena, =0, equation (11) describes a nonlinear oscillator with a nonlinear positive 
damping and a nonlinear restoring term, both of them without time lags. In this case the 
steady state is globally and asymptotically stable and sustained oscillations of thermal power 
are impossible. Ifa, #Oan additional restoring term appears, with a continuous distribution 
of mechanical time lags given by the memory function K (1 ), so that the possibility of 
sustained power oscillations can not be excluded on a priori grounds. 


2.3. Approximation by a Discrete Delay: Power Thresholds and Fast Instabilities 


Following a well grounded tradition in control theory and electrical engineering, in which the 
emphasis is overwhelmingly on discrete lags, let us try to approximate the continuous 
distribution of lag times in equation (11) by a single discrete lagr,: 


+ foe? 9(x)+ 0%, -g(alt-1,))=0 (14) 
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The discrete lag should be chosen such that the qualitative behavior of equations (11) and 
(14) is essentially the same in the different scenarios of interest. This may be accomplished 
through an average lag ¢ and a tuning dimensionless factor 0 , which is of the numerical order 
of 1. 

The average time lag may be estimated using the non-negative damped exponential: 


— Sin Ont 
This exponential describes the mechanical damping of the structural materials of the reactor, 
for sub-critical, just-critical and super-critical damping. The impulsive response function K (: ) 


is not adequate for time averaging purposes, because if the system is under-damped it 
oscillates changing its sign. 


oO 


The normalizing factor ¢,,,-@,, 18 introduced in order to have J E (¢ ). dt=1 


1, = 0-7 =0-[t-E)-dt =8. =0- (15) 
0 


We see that the time lag increases with decreasing mechanical damping. 
A linear approximation to (14) in a neighborhood of the steady state (x =0, < =()), taking 
t 


into account equations (5) and (6) gives: 


eG. + 6p) MO oy? xl) +0%, 20 (16) 


The stability of the steady state can be studied, as usual for linear delay differential equations, 
using the ansatz 


x(t) = e*" (17) 
Substituting this tentative solution in (16), it follows a transcendental equation for z : 
2 2 2 —2tg 
Zz +(cy te, )z+o, +a@°y-e “4 =0 (18) 


As the left hand member of (18) is always positive for z positive or zero, there are no positive 
or zero roots. If the system gets unstable, it must do so at threshold oscillating with a certain 


frequency @, such that z= j-q@,is a root of (18). (Here j = /—1). Taking the steady state 
power P, as control parameter, an analysis of the roots near the critical value of P, and the 


stability of the steady state at the critical value show that this corresponds to a supercritical 
Hopf bifurcation ([7]) producing a stable limit cycle when the steady state gets unstable for a 
critical value of the delayt,. For z-t,small enough:e““ ~1—z-t, Introducing the 


combined nuclear frequency @, , equation (18) may be approximated thus: 
2 2 2 
22 +(cy +c, —@’y -ty )z+a'n XO (19) 


Oy = Or a @’y (20) 
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From (19) together with (6a), (6b) and (12) it follows that an approximation to the critical 
mechanical delayr, . and to the frequency of oscillation w, at the instability threshold are: 


: winter BC 1, [ail (21) 
ae wy ayb Po ayb 
QO, © Oy (22) 


o 
“is small enough, we can expect that the linear approximation to 
m 


2 
So, if Oy ‘tg =0.- 


ee , may be good enough for our purpose. 
If the delay is greater thant, .we should have power oscillations around an unstable steady 
State. 


a,|C 
The critical delay at very large powers reaches its minimum value —— . 
a 
y 
Then, if the mechanical delayt, is less that this minimum, the reactor will be always stable 
independently of its steady state power. 


When |@,|= Obuta, #0 (as in Thompson’s model), equation (21) shows that no matter how 


small is the mechanical delay t,, there will be a steady state power above which the reactor 


Qa; 


will always be unstable. 
Now, givent,, from equation (21) we obtain an estimation of a threshold P, to power 
instability, valid if the prompt feedback coefficient 
B-C 
P,= eae (23) 
i- a,|°C 


a, -b-t, 


ar,| is small enough: 


Ift, is near enough tot, ., the period of the oscillations could be estimated by 
20 


2 
Cy +c ao ft 
y 
Ow N F ? d 


(24) 


3. CONCLUSIONS 


- (1) In the extension of the model of Weinberg-Wigner-Thompson proposed in reference [1] 
and discussed in the present paper, a nonlinear nuclear-thermal oscillator (that corresponds to 
reactor point kinetics with thermal-elastic feedback and with frozen delayed neutron effects) 
is coupled with a linear mechanical-thermal oscillator (that corresponds to the first normal 
mode of mechanical vibrations excited by thermo-elastic effects). The nonlinear damping of 
the nuclear-thermal oscillator, given by (5a), is always positive (stabilizing). In reference [2] 
the results of several digital simulation runs are given for a model without prompt power 
feedback. Analogous runs should be done with the model proposed in the present paper. 
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Chaotic behaviour could be produced for small enough values of the prompt feedback 
coefficient of reactivity. 


- (2) A nonlinear integro-differential equation was derived to describe thermal power kinetics, 
assuming that the mechanical oscillator was at rest when a sudden perturbation in the reactor 
power was produced. It could be interesting to study what happens with the reactor when a 
sudden perturbation is introduced in the mechanical oscillator while its power is in its steady 
state value. 


— (3) The continuous distribution of mechanical time lags was substituted by an equivalent 
discrete time lag. The stability of the linear approximation to the corresponding nonlinear 
second order delay-differential equation was studied in a neighbourhood of the steady sate of 
the reactor. 

In reference [1] a modification of a first order delay differential equation for thermal power, 
proposed by Ackasu, Lellouche and Shotkin ([8]), was posed to describe fast reactor kinetics 
with prompt and delayed power feedbacks. There the discrete time lag was introduced as an 
independent parameter. Here the discrete time lag is obtained as function of the parameters of 
the mechanical oscillator and is a delay exclusively of mechanical origin. 


- (4) When the steady power verifies P, > P,, with P,, given by equation (23), the coupling 
between the nuclear oscillator and the mechanical oscillator produces fast (relative to 
delayed neutrons time scale) and sustained power oscillations around an unstable steady 
State. 


- (5) The threshold of power P,, doesn’t exist if the prompt power feedback is strong enough 


Dy 


a,-b-t 
anda, verifies the inequality |a,| > that stems from (23). According to this model, 


a; 


dangerous power oscillations may appear only if the inequality is reversed. 


- (6) These results could be of some interest in order to discuss with engineering educated 
people that opposes to the expansion of nuclear power technology, arguing with some aspects 
of the physical safety of nuclear reactors that are indeed outside the scope of most available 
neutronic-thermal-hydraulic numerical codes. 


- (7) The properties of the point kinetics model with feedback proposed in this paper (and in 
[1]) could be studied further from the perspective afforded by the powerful and well 
developed asymptotic methods for the analysis and control of nonlinear systems [9]. 


- (8) Point kinetic models with constant reactivity feedback coefficients, like the present one, 
are often too crude idealizations for studies in reactor dynamics and control, so that detailed 
numerical and experimental research for specific reactor types must be done to discover the 
many subtleties that seem to be hidden under the umbrella of nuclear-mechanical coupling. 


- (9) An intermediate step, between a somewhat oversimplified point kinetics modeling of 
nuclear mechanical coupling and detailed numerical or experimental studies with specific 
reactor types, could be a still mainly analytical but distributed parameters approach to reactor 
and NPP stability and control like the one that appears in references [10] and [11]. 
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An analytical approach to bifurcations and stability in 
simplified mathematical models of nuclear reactors 


Roberto Sudrez-Antola 


ABSTRACT 


Asymptotic analytical methods are applied to study some problems related with bifurcations (both local and 
global) and stability in three simple mathematical models of nuclear reactors. 

The first case is a reactor which is subcritical at rest and driven by a distributed external neutron source. Under 
suitable conditions in the temperature feedback reactivity, the mathematical model predicts a static global 
bifurcation with three critical states, two stable and one unstable, with the possibility of a runaway. The dynamics 
of this system is studied, in a framework of slow manifold theory, by methods of restricted nonlinear modal 
analysis, and the results of a digital simulation are summarized. The model could be modified and extended to 
study the space time dynamics of sub-critical multiplying systems driven by external neutron sources (neutron 
beams produced by accelerators). 

The second case is related with in phase and out of phase xenon oscillations in large thermal reactors. The known 
subcritical Hopf bifurcation that appears in the context of global mode oscillations is revisited. After applying a 
nonlinear modal analysis to the mathematical model, a normal form is derived by an averaging method. 
Approximate analytical formulae for the radii of the unstable limit cycles and for the trajectories of the state 
variables in a neighborhood of the bifurcation point are obtained. 

The third case is a non-trivial modification and development of a simple mathematical model intended to describe 
certain mechanical kinetic effects stemming from a possible coupling of nuclear, thermal and mechanical vibration 
processes. We show that, under suitable conditions, a dynamic supercritical Hopf bifurcation in the reactor power 
appears in the framework of our modified model. A normal form is derived by an averaging method. Analytical 
formulae for the radii of the stable limit cycles and the trajectories of the state variables in a neighborhood of the 
bifurcation are given. 


Keywords: nuclear reactor dynamics, reduced order models, restricted nonlinear modal analysis, 
averaging methods, bifurcation, stability, mathematical modeling, homogenized equivalent reactor 


1. Introduction 


One of the goals of reactor design and operation is to restrict the possible states that the reactor can 
reach, during steady operation and during transients. For example, after a disturbance that occurred 
during the steady state operation of a reactor, it is desirable that the perturbed state of the system does 
not get too far away and returns quickly enough to its original steady state. Also, during transients, 
certain restrictions must be imposed on the time scales of change in the reactor's state. 

The modification of certain parameters of the reactor (like steady power or coolant flow) can destabilize 
the steady state. To avoid it is necessary to identify the instability thresholds in the parameters of the 
reactor, beyond which the steady state loses its local stability and its dynamic behavior suffers a 
qualitative change, like the neutron-thermal-hydraulic oscillations that under certain conditions appear 
in boiling water reactors (BWR) (Turso et al., 1997; Rizwan-uddin, 2006), or the xenon oscillations that 
can appear in large reactors (Henry, 1975; Duderstadt and Hamilton, 1976; Lewins, 1978) 

But even if the steady state is stable under small perturbations, it can be unstable when the amplitude 
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of the perturbation is large enough. Once this threshold amplitude is exceeded, the state of the reactor 
will move farther away, perhaps approaching to a new and undesired steady state or being attracted 
towards an oscillating pattern in the state variables. 


To study in detail these stability problems in nuclear power plants (NPP), in large radioisotope 
production reactors or in research reactors, full scale realistic mathematical modeling with focus in the 
reactor must be done beginning with a suitable coupled system of non-linear partial differential 
equations. Then, the corresponding computational code for digital simulation of steady states and 
transients is constructed and tested. At this level of analysis random fluctuations in reactor parameters 
are regarded as a nuisance that contaminates the deterministic evolution of the state variables, like 
radiation background in detectors, and are neglected. At a deeper level of analysis, the fluctuations must 
be included in the mathematical model, because under certain conditions, random fluctuations can shift 
the state of the reactor out of its steady state stability region and modify the stability thresholds 
mentioned above, amongst other effects that are outside the scope of deterministic mathematical 
models. (Konno et al., 1992; Konno et al., 1994; Hennig et al. 2019 (a), this issue). 


The objectives of the rest of this introduction are: 

1.1- To review some basic topics of dynamic systems, both deterministic and random, needed as a 
framework for the present research. 

1.2- To introduce three asymptotic methods (inertial manifolds, central manifolds and slow manifolds) 
that can be used to reduce the number of state variables in mathematical models of nuclear reactors. 
1.3- Establish the specific research objectives of this paper. 


1.1- Review of some basic topics of dynamic systems 


1.1.1-Deterministic dynamic systems 

Considered as a deterministic dynamical system, the above-mentioned dynamic equations of the 
reactor, that include partial differential equations, describe the movement of a point x (the state of the 
reactor) in an abstract space_X of infinite dimensions, the phase space of the system. This movement 
follows an evolution equation that can be written like this (Attle-Jackson,1989): 


< - F(x) 


Here F represents a nonlinear operator that includes partial derivatives relative to the Euclidean spatial 
coordinates, E=(G5enye;) is a vector whose components are fixed reactor parameters and ft 
represents, as usual, the time instant. In the so-called reduced order models (ROMs) of nuclear 


reactor dynamics, the space_X is of finite dimension and the evolution equation is a set of nonlinear 
ordinary differential equations. 


Each set of possible values of the parameters c = (c, pCpsecs Cy ) can be considered as a point belonging 
to a k-dimension Euclidean space C= R*: c= (esGs kee ) eR Wecall C=R* the parameters 
space. 

Given enough regularity, for each initial condition x, there is a unique trajectory x(t »Xo3 c) . The set 


of phase points that belongs to a trajectory determine an oriented curve in phase space name orbit. The 
family of all these curves is the phase portrait of the dynamical system for the chosen set C,,C,,...,C; 


of parameter values. 
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When a change in some parameter or parameters (called control parameters) produces a qualitative 
(topological) change in the phase portrait, it is said that a bifurcation has occurred. 


For example, a stable equilibrium (fixed) point ra (ca) is destabilized and a stable limit cycle of 


oscillation appears (a kind of local bifurcation of the dynamic type) or new equilibrium (fixed) points 
appear at a certain distance from the already existent equilibrium points in phase space without changing 
the type of pre-existing fixed points (a kind of global bifurcation of the static type). 


The bifurcation diagram of a deterministic dynamical system is the partitioning of the parameters 
space into regions where the phase portraits are qualitatively (topologically) equivalent. 

A set included in phase space, such that if a point of a trajectory belongs to the set, then all the points 
of the trajectory belong to this set, is called an invariant set. (A fixed point, the points of a cyclic 
trajectory and in general the points belonging to a set of trajectories is an invariant set). 

Sometimes invariant sets, regular enough to be considered as a refinement and a generalization of the 
regular curves and surface of ordinary Euclidean space EF, can be found in the phase space of a dynamic 
system. This kind of sets, known as finite dimension invariant manifolds locally resemble an Euclidean 
space E, (the same E,, everywhere, son is the dimension of the manifold), but globally it might be 
more complicated. (For example, a torus and a Klein bottle are two-dimension manifolds). 

An attractor is an invariant set such that all trajectories passing through nearby points tend to that set. 
(Stable fixed points and stable limit cycles are attractors). 


As will be seen below in relation with the mathematical foundation of the reduced order models, certain 
higher dimension invariant manifolds can be attractors. 


1.1.2-Random dynamic systems 

Although the theory of dynamic systems and one of its branches, bifurcation theory, ignore fluctuations, 
as Hermann Haken (Haken, 2003) has repeatedly emphasized fluctuations are crucial precisely at and 
in the close neighborhood of the bifurcation thresholds. If random fluctuations in the parameters are 
incorporated in the mathematical model, a random dynamical system is obtained, described by a 
random differential equation: 


5 wf (Ber 8e() 


Here 6¢ (t) represents the random fluctuations of the reactor parameters, X(t) is a stochastic process 


: dx. . : Seog 
that represents the random state of the reactor at each instant of time and ae is its stochastic derivative. 
t 


(van Kampen, 2001) (Gardiner, 2004) (Scott, 2013). 
In the limit of large numbers of particles involved in random events in a representative volume element 


of a reactor, the random state can be expressed as the sum of a deterministic component x(t) and a 
random component OX(t): X(t) =x(t)+ X(t) 
This approach is known as the linear noise approximation (van Kampen, 2001; Scott, 2013). 


In the reactor case we assume that the term x(t ) is a not chaotic solution of the usual full scale 
deterministic equations of nuclear reactor dynamics, and the fluctuating term OX(t ) is described by a 
probability density distribution ITI (Ox, t) centered upon x(t) . It is possible to imagine that II (Ox, t) 


moves along the deterministic term x (t) . 
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When, as in the case of reduced-order models, the phase space is of finite dimension n , the fluctuation 


term is composed byn scalar random processes 5%, (t),5%, (t),...,0%, (¢), the corresponding 
multivariate probability density 11 (6%, (t VOX (t ), si OX, (t);1) is Gaussian and is obtained solving 


a Fokker-Planck equation that describes the evolution of the probability density of the fluctuations. 
(Gardiner, 2004; Scott, 2013). 


Away from the bifurcation points the width of the probability distributions of the state variables are in 
general small relative to their mean values. This allows the estimation of the deterministic values of the 
parameters of the reactor, from measured data, by well-known statistical methods (Gershenfeld, 1999). 
When the linear noise approximation can’t be applied, there are more powerful methods that in principle 
allow us to distinguish between chaos, quasi-periodicity and random noise, as well as to identify 
possible attractors of the dynamic system. These methods are based on embedding techniques and the 
construction of information dimensions (Suzudo et al.,1993; Gershenfeld,1999). 


As stressed by several authors (Konno et al., 1992; Konno et al., van Kampen, 2001; Gardiner, 2004; 
Scott, 2013) random fluctuations can modify the bifurcation scenario. 

However, under conditions of applicability of the linear noise approximation, a deterministic analysis 
can explain the stability behavior. An example are the bifurcations found for the NPP Leibstadt (Hennig 
et al., 2019 (a), this issue). 


1.2- Reduction of the number of state variables in mathematical models of nuclear reactors: 
methods of ROMs development. 


Let us return now to the deterministic approach to nuclear reactor dynamics. The non-linear partial 
differential equations that appear in full scale deterministic mathematical models of nuclear reactor 
dynamics belong to the class known as dissipative nonlinear partial differential equations (DNPDEs). 
This kind of equation describe mechanism of energy relocation and dissipation in the reactor, as well 
as the interaction between these mechanisms, leading to the appearance of complex dynamics such as 
the nonlinear oscillations mentioned previously. 

Although there are theoretical methods that allows to obtain certain general quantities related with the 
solutions of DNPDEs and its asymptotic evolution in their corresponding phase spaces (Logan, 1994; 
Chueshov, 2002), these results, very interesting from a mathematical point of view, often don’t give the 
details usually requested by physicist and engineers. 

However, many theoretical and computational studies have shown remarkable similarities between the 
long-time (asymptotic) evolution of the solutions of certain dissipative nonlinear partial differential 
equations (DNPDEs) in their corresponding infinite dimension phase spaces and the solutions of certain 
ordinary non-linear differential equations (ODEs) in their own finite dimension phase spaces. 
(Malinietski, 2005; Hennig et al. 2019(a), Hennig et al. 2019 (b), this issue) 

Furthermore, the approach of the solutions of the DNPDs to the solutions of the corresponding system 
of ODEs is in many cases fast enough to allow the dynamic study to be done entirely using these ODEs, 
reducing the computational effort. (Constantin et al. 1989; Malinietski, 2005) 

Basically, the mathematical background that justifies the introduction of reduced order models (ROMs) 
in studies of NPP dynamics and control is the existence of hidden attractors in phase space (a stable 
fixed point, a stable limit cycle or a more complex finite dimension attracting manifold). 
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1.2.1-Method of ROMs development by inertial manifolds 

The origin of this similarity in the asymptotic behavior of DPDEs and ODEs is due to the following 
fact: the asymptotic evolution of the solutions of these DPDEs occurs near and exponentially 
approaching to a finite dimensional invariant attracting manifold in phase space, the so-called inertial 
manifold. For each trajectory that begins in a neighborhood of the inertial manifold, there is a trajectory 
included in the inertial manifold such that the distance between these two trajectories tends to zero 
exponentially. These and other properties of the inertial manifolds are described by several authors 
(Constantin et al. 1989; Chueshov, 2002; Malinietski, 2005). 

Now, this inertial manifold corresponds to or is very well approximated by the solutions of a suitable 
set of nonlinear ordinary differential equations: with the inertial manifold theory, it is shown that certain 
dissipative PDEs have the same asymptotic behavior of a finite dimensional ODEs system, known as 
the inertial form corresponding to the considered inertial manifold. 

Furthermore, the exponential approach of the solutions of the DNPDs to the corresponding inertial 
manifold (solutions of the ODEs) is in many cases fast enough to allow the dynamic study to be done 
entirely on the inertial manifold, that is, using ODEs and reducing computational effort. As 
consequence, if an inertial manifold exists, in principle at least, a reduction of the dimension of the state 
vector of a complex nuclear system can be done, preserving the main properties of the original model 
and attaining a small approximation error in the framework of the simplified mathematical model. 
(Hennig et al. 2019 (a) and 2019 (b), this issue). 

Unlike other asymptotic methods, in principle at least, inertial manifold theory could be applied to 
simplify the analysis of a dynamic system at all points of its parameter space: a reduction by inertial 
manifold (IMR). 

However, the construction of the exact inertial forms that correspond to a given system of DNPDs is 
often not feasible. So, in practice only approximate systems of ODEs (approximate inertial forms) can 
be constructed to reduce the dimension of the phase space of a complex system (Hennig et al. 2019 (a), 
this issue; Manthey et al. 2019, this issue). 


Besides inertial manifold theory (a global method), there are two general procedures to reduce a high 
dimension dynamic to a low dimension dynamic: the center manifold theory (a local method) and the 
slow manifold theory (a global method) (Wasow, 1965; Nicolis, 1995). 

In both cases a few dominant mode amplitudes can often be identified, so that the behavior of the system 
is essentially determined by these dominant modes. 


1.2.2- Method of ROMs development by center manifolds deterministic and stochastic 
In the center manifold theory, the stability of the equilibrium solutions (fixed points) x (c) in phase 


space X is studied in parameters space C = R* varying systematically some parameter values (control 
parameters). 
The system is expressed in each neighborhood of an equilibrium solution as the sum of two terms, one 


linear J(X(c))-(x-x¥(c)) and the other non-linear h(x(c);x-x(c)). This last one vanishes 
faster (at least quadratically) than the linear term when the state of the system x tends to the equilibrium 


point x(¢} in phase space. So, we have: 


& — 3(x(c))(x-(6))+4(F(€):x-2(6) 


For a nonlinear system of ordinary differential equations of order n, X = R”. 


Original version: Progress in Nuclear Energy 114 (2019) 171-190 Special Issue on Nonlinear Stability 
Analysis: Nonlinear nuclear reactor stability analysis and special aspects of modelling of complex 
dynamical systems. 


The corresponding fixed points X(c) are called hyperbolic (or non-degenerate) when all the 


eigenvalues of their Jacobian matrix J (x (c)) have non-zero real parts. This occurs almost always in 


parameters space. 
According to the Hartman-Grobman theorem, also called principle of linearized stability 
(Guckenheimer and Holmes, 1983; Attle-Jackson,1989; Nicolis, 1995; Hennig et al. 2019 (a), this 


issue), the behavior of the dynamic system in a neighborhood of a hyperbolic point x (c)is 
; Caer dace ot ee a = 
topologically equivalent to the behavior of the linearized approximation a =J (x (c)) : (x —xX (c)) : 
t 


However, for certain points c,,, (bifurcation points or critical points) in parameters space C=R* the 
corresponding equilibrium point x (Gy) in phase space X =R”" can be non-hyperbolic (or 


degenerate): they have some eigenvalues with zero real part. In case of complex conjugate eigenvalues, 
there is at least a pair of complex conjugate eigenvalues laying on the imaginary axis of the complex 
plane. So, the bifurcation points (critical points) in parameters space, where the real part of one or in 


general a few eigenvalues of the linearized system in a neighborhood of x (c) change its sign, are 


identified. 
dx = = = = : : 
In the local development ah = J(X (c)) . (x —X (c)) +h (x (c) 5X-X (c)) of the dynamic equations, 


the nonlinear term h(x (c);x —XxX (c)) can't be left aside when one seeks to build a dynamic system 


topologically equivalent but easier to analyze than the original. 


In a neighborhood of these critical values of the parameters, where the equilibrium points lose their 
stability, the change of stability is often preceded by the critical slowing down of a few mode amplitudes 
(the dominant ones) that slave the others in such a way that (perhaps after a short transient) the evolution 
of the state of the system occurs in a low dimension manifold (which exactly at the critical point is the 
so called center manifold in bifurcation theory) (Guckenheimer and Holmes, 1983; Nicolis, 1995; 
Haken, 2003). 

The other modes of evolution are given as functions of the dominant ones, called order parameters in 
Synergetics (Haken, 2004). 

This reduction in the description of the dynamic system to a dynamic in terms of the dominant modes 
of evolution is known as center manifold reduction (CMR) and generates a reduced order model (ROM) 
linked to the CMR. 

By near identity transformations of state variables (Guckenheimer and Holmes, 1983; Nicolis, 1995; 
Kuznetsov, 1998), in a neighborhood of the bifurcation point it is possible (although not always easy) 
to build a dynamic system that gives a local description of the dynamics topologically equivalent to the 
original system, but easier to study. The nonlinear differential equations of these topologically 
equivalent dynamic systems are known as normal forms and allow a classification of the local 
bifurcations. 

In principle the normal forms related with CMR are asymptotically exact, unlike the inertial forms 
related with RIM accessible in practice, which are asymptotically approximate. 

The concept of a normal form as an asymptotic local description of the bifurcation dynamics of a system 
in the neighborhood of an equilibrium point, is not exclusively related to the classical approach initiated 
by Poincare, the method of multiple time scales or the more recent variant developed by Arneodo, 
Coullet and others (Nicolis, 1995). Normal forms can be obtained also by the averaging methods 
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described below, in subsection 2.4. 


Local bifurcation theory has several advantages to cope with stability problems in nuclear reactors: 
many stability issues can be understood as a change in the stability condition of a steady state operating 
point, including the birth of limit cycles of oscillation (Hopf bifurcations) either stable (supercritical) 
or unstable (subcritical). 

But other stability problems seem to be outside the span of local bifurcation theory. For example, the 
derivation of the boundaries of instability around a stable steady state (operating point) in case of having 
at the same time at least another stable steady state with the possibility of an undesired run away from 
the operating point to this other steady state. A Bautin bifurcation of limit cycles is another example, 
well studied recently in relation with BWR stability problems (Lange, 2009; Lange et al, 2011; Lange 
et al, 2013). 


When random fluctuations can’t be neglected, it is possible to resort to a random center manifold 
theory, which allows a stochastic approach to bifurcation theory (Konno et al., 1994). 
Random normal forms appear, in general as nonlinear Langevin equations like this one, in terms of the 


random distance a (t) from a steady state and a random phase (t) : 


—a(t)=K-a(t)-y-@ (t)+N, (t) <#(t)=0(a)+8, (1) 


Here N, (t) and N, (t) are random noises, «is a real and y is real and positive. 

Without the noise, the corresponding deterministic normal form describes a supercritical Hopf 
bifurcation that occurs when « changes its sign from negative to positive. In presence of the noise this 
formalism describes a stochastic supercritical Hopf bifurcation. 

If the linear noise approximation is valid, usual statistical methods can be applied to estimate the 
numerical values of parameters, in this last case « and y, from noisy data, as mentioned at the end of 
1.1.2. 


In a neighborhood of the bifurcation point, when linear noise approximation usually fails, but N a (t) is 


a Gaussian white noise of zero mean and known variance o”, it is possible to derive a formula 


for the probability density function (pdf) of @ (t ) averaged over on cycle of oscillation (Konno et al., 
1992). The only parameters that appear in this formula are - and a (unless numerical factors), 
Oo Oo 


so if o is known, the other parameters could be estimated from averaged random data. 
Besides, it is possible to apply the more powerful methods framed in embedding techniques 
and information measures mentioned at the end of 1.1.2. 


1.2.3- Method of ROMs development by slow manifolds 

The slow manifold theory, when applicable, allows us to identify a few dominant degrees of freedom 
in phase space, even when the equilibrium solutions are far from bifurcation points in parameters 
space. A so-called slow manifold, of low dimension, has to be identified in phase space. Furthermore, 
this slow manifold must be stable in the sense that any trajectory located in a neighborhood of the 
manifold, after a short transient, should approach to and remain near this manifold. It is called slow 
because after this approach to the manifold, the state of the system changes slowly in comparison with 
the initial transient. 
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According to this definition, the center manifold defined in 1.2.2, on which the relevant part of the 
dynamics lies, is an example of slow manifold, because fast modes of evolution are slaved by slow 
modes in a neighborhood of a bifurcation point. But not every slow manifold is a center manifold (see, 
for example, Nicolis, 1995): it is not necessary to have a bifurcation point of the dynamic system in the 
set of parameters values considered in slow manifold analysis. 


Moreover, a slow manifold is not necessarily formed by orbits of the dynamic system (a slow manifold 
is not necessarily an invariant manifold like the inertial manifold is), although some orbits must pass 
close enough to the slow manifold for a certain interval of time at least (Wasow, 1965; Tihonov et al., 
1970). 


Now we are going to summarize some relevant aspects of Tihonov’s approach to singular perturbation 
theory, which gives a foundation to identify and characterize slow manifolds. 

Let us suppose that the equations of evolution of a deterministic dynamic system can be written as 
follows (Tihonov et al., 1970): 


eS =F (z,y0) “— F(z,y:0) 


The state variable is given now by the pair of vectors (2 y) interconnected through the preceding 
evolution equations. As before, c represents a point in parameters space and the new positive parameter 
€ is small enough and in the limit will equal zero. This means that z(t sc} varies much faster than 
y(tsc) when 2 is not too close to F (z, y3c) = (0). The initial conditions of the given trajectory in 
phase space are 2(tj3¢) =Z, and y(t);c) =V- 

Now we consider the degenerate dynamic system whose evolution equations are obtained from the 
original ones putting ¢=0: 

F(z, y;c)=0 —= f(z, yc) 

In general, the solution of the equation F (z, y3c) =QOis a set of branches 
z=h, ( y3c) (A =1,2,...,m) that do not intersect and can be classified as stable branches 


z=h,(y;c) and unstable branches z=h,,(y;c) A branch h(y;c) is stable if all the eigenvalues of 
OF . . 
the operator Fe il y;c) ; y;c) have negative real parts (Tihonov et al., 1970). 
Zz 


These stable branches are the stable slow manifolds mentioned above: given initial conditions belonging 
to the region of attraction of the stable branch, the trajectories z (¢ 5 c) and y (¢; c) first approach to 
the stable branch and then remain near it. 

Since z(t;c) changes much faster than y(tsc) as they approach the stable branch, it is possible to 


introduce the following simplification of the mathematical model, known as the inner approximation, 
keeping the slow variable fixed at its initial value: 


d 
e— =F (2,956) 


This simplified equation is solved with the initial condition Z (3c) =Z, and an inner solution 
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ae (¢; c) can be constructed. 


Once near the stable branch, a second simplification of the mathematical model is obtained, known as 


the outer approximation, substituting z= h, ( y; c) in the equation of the slow variable: 
dy 
—= f(h,(y:¢),¥¢) 
dt 
This new simplified equation is solved with an initial condition y, that verifies z, = h, ( Vis c) , and an 


outer solutions Vig. (t3C) ANd Z yu, (56) =A (Vouer (6 c);c) 


The inner and outer solutions can be combined to obtain a valid approximation including the fast 
relaxation from the initial conditions to a neighborhood of the stable manifold, and the slow movement 
near the manifold (Tihonov, 1970; Lin and Segel, 1988). 

The slow manifold theory summarized above may be applied to study transients, including nonlinear 
oscillations, in nuclear reactor dynamics, because in a complex nuclear system there is a hierarchy of 
widely separated time scales. 

This goes from prompt neutron — dominated effects (hundredths of seconds), to heat transfer from fuel 
to coolant (tenths of seconds), coolant transit times through core and precursor dominated effects 
(seconds), coolant transit time through the entire primary circuit (tens of seconds), diurnal electric load 
variations and xenon flux tilting (tens of hours), samarium production (months), and fuel burn-up and 
transuranic isotope production (years). 

When the variables of interest (variables of reference) belong to a certain time scale, and there is an 
attracting slow manifold, it is possible to simplify the dynamic analysis applying the following two 
principles to link scales of different orders of magnitude: 

1-The variables belonging to processes with time scales at least an order of magnitude greater than the 
reference time scale, can be considered as frozen. This corresponds to setting the values of the slow 
variables to build an (approximate) inner solution (in Tihonov’s sense) for the fast variables. 

2- The variables (fast) belonging to processes that evolve with time scales at least an order of magnitude 
smaller than the variables of reference (slow), after a short transient (produced in the so-called inner 
time scale) can be considered as relaxed to equilibrium with the reference variables (evolving in the so- 
called outer time scale). This corresponds to the building of an (approximate) outer solution (in 
Tihonov’s sense) for the slow variables. 

Using these principles, the slow manifold and a few dominant degrees of freedom of evolution, each 
with its corresponding amplitude (a function of time), can usually be identified. The dominant degrees 
of freedom evolve more slowly and tending to slave the evolution of the other degrees of freedom. So, 
an infinite dimension dynamic or at least a finite dimension dynamic with a big number of dimensions, 
is reduced to a low dimensional one: a ROM is obtained. 

Often slow manifold theory can be applied far from bifurcation conditions, to further reduce the 
dimensionality. 


1.2.4-ROMs in nuclear reactor dynamics 

The use of ROMs in nonlinear stability analysis of nuclear reactors allows the application of a suitable 
combination of semi-analytical bifurcation theory (including software specific to study bifurcations) 
with digital simulation of the dynamics. The results thus obtained can be used to guide the numerical 
simulation using full system codes for DNPDE. 


A typical ROM for a BWR may have a number of state variables near 20, interrelated by a 
corresponding system of near 20 nonlinear ordinary differential equations. 
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There are simpler mathematical models, like the classical March-Leuba’s model, that have few state 
variables (5 state variables in March-Leuba’s model) (March- Leuba et al., 1986). 

However, these highly simplified mathematical models reproduce some of the main characteristic of 
nuclear reactor’s dynamics, are easier to study and can be derived, making suitable simplifications, from 
a more comprehensive ROM. So, this kind of simple model could be considered a ROM of a ROM. 
When the ROM is simple enough, it is sometimes possible to study local bifurcations of steady states 
and limit cycles, and even some global bifurcations, using approximate analytical methods. The results 
thus obtained usually offer a useful guide to deeper stability and bifurcation studies using the same 
model or more complex ROM’s of the same physical system. 


1.3- Specific research objectives 


The purpose of the present paper is to apply asymptotic analytical methods to study some problems 
related with bifurcations (both local and global) and stability in three simple ROMs of nuclear reactors. 


(a)-The first case is a reactor which is subcritical at rest and driven by a distributed external neutron 
source. Under suitable conditions in the temperature feedback reactivity, the mathematical model 
predicts a static global bifurcation with three critical states, two stable and one unstable, with the 
possibility of a runaway. The dynamics of this system is studied by the method of restricted nonlinear 
modal analysis. (This method is introduced and explained in subsection 2.3. below). A saddle-node 
bifurcation is found in the context of a static global bifurcation. Some related digital simulation results 
are summarized. 


(b)-The second case is related with in phase and out of phase xenon oscillations in large thermal 
reactors. The method of restricted nonlinear modal analysis is applied to a simplified model of the 
reactor. The known subcritical Hopf bifurcation that appears in the context of global mode oscillations 
is revisited. A normal form for this dynamic bifurcation is derived by an averaging method introduced 
and explained in 2.1.3 below. Approximate analytical formulae for the radii of the unstable limit cycles 
and for the trajectories of the state variables in a neighborhood of the bifurcation point are obtained. 


(c)-The third case is a non-trivial modification and development of a simple mathematical model that 
appears in a very recent book of Physics of Nuclear Reactors (Marguet, 2017). This model is intended 
to describe certain mechanical kinetic effects stemming from a possible coupling of nuclear, thermal 
and mechanical vibration processes. We show that, under suitable conditions, a dynamic supercritical 
Hopf bifurcation in the reactor power appears in the framework of our modified model. A normal form 
is derived by the averaging method mentioned above and described in 2.1.3 below. Analytical formulae 
for the radii of the stable limit cycles and the trajectories of the sate variables in a neighborhood of the 
bifurcation are given. 


Now, after having specified our research objectives, in the following section we review some of the 
mathematical tools necessary to achieve these objectives. 


2. Methods 


The three simple ROMs mentioned above can be derived from the general mathematical model of an 
equivalent homogeneous reactor summarized in section 2.1, after making one of the approximations 
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described in section 2.2. The model incorporates mechanical variables related with the structural 
elements in the core, needed to build the third ROM. 

A variant of the method of nonlinear modal analysis, needed to develop the first two ROMs and to 
analyze the global static bifurcation in the first ROM, is summarized in section 2.3. 

An averaging method for nonlinear differential-integral equations, needed in the analysis of limit cycle 
bifurcations in the last two ROMs is reviewed in section 2.4 


2.1. General Background 


A fairly general, albeit already simplified mathematical framework that can be used to describe the 
dynamics of a nuclear reactor, including a power reactor belonging to a NPP, is given by the following 
equations, jointly with suitable initial and boundary conditions: 


10 . . ‘ 
+P = (1-p)M[d.xy,w]-L[bx.y.n] + Ac +8 ‘ 
k=1 

a) - 

= hve [9 eM w] AC, 0) 
a : ~ ° 4 (3) 
> 7 vy, Wy w,, és 
= = w\¢, Ww, Ww, Z 
Ow, = W,, fw, w, | ‘ 
Ot 


As usual, time instants will be represented by ¢ and spatial points in the core by position vectors 7 . 
The first three equations describe the neutron kinetics in the core, using a one group diffusion 
approximation, being ot, 7) is the neutron flux. A multi-group description is not considered here, to 


allow a simpler analysis. 
An external distribution of neutron sources S (¢, 7 ) appears as the last term in equation (1), which is the 
balance equation of the space-time neutron distribution. 


The fields c, (t, r Jare the concentration of delayed neutron emitters. Their kinetics is represented by 


the equations (2). As usual, /, represents the fraction of the k group of delayed neutrons, 2 = >» P, 
k 


the total fraction of delayed neutrons and 4, the decaying constant of the k group of neutron emitters. 


The scalar fission operator M and the scalar absorption and leakage operator L are linear operators 
acting on the field dt, Fr i In general, these operators will be functions of the feedback variables x , y 
andw. 

Because the approach to reactor dynamics in this paper is analytical, the equivalent homogeneous 
reactor corresponding to the real non-homogeneous one will be considered from now on. So, we work 
with a fictitious mix of fuel and moderator with effective macroscopic cross section of fission, 
absorption, scattering and transport at each point of the core (Stacey, 2018). 
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In the one group diffusion approximation in the equivalent homogenized reactor adopted here, if v is 


the average number of neutron produced per fission and X , is the effective macroscopic fission cross 
section, the production operator is given by the formula: 
M[g|=v 2-9 7) 
The loss operator is given by the sum of a neutron absorption term with effective macroscopic 
absorption cross section 2, and a local leakage term given by the divergence of the neutron current J 
(in this case e represents a scalar product): 
Lg|==,-¢+VeJ (8) 

If D is the effective neutron diffusion coefficient and V @ is the gradient of the neutron flux, the neutron 
current is given by: 

J=-D-V¢ (9) 
Then: 

[|¢|=2, -¢-Ve(D-V¢) (10) 

As the diffusion coefficient of fuel is not very different from the diffusion coefficient of moderator, and 
the volume of moderator is usually fairly large compared to the fuel volume, an effective diffusion 


coefficient equal to the diffusion coefficient in the moderator can be introduced to simplify the 
mathematical model of a heterogeneous thermal reactor. If, in a first approximation, the diffusion 


coefficient is considered space-independent, the loss operator simplifies to L\¢|= x, :9-D: Vo 


and equation (1) can be further simplified and re-written as follows: 


1 og 24 
Tap ABV Ey - 2a Gt DV i eae (11) 


The effective parameters 2 ,, 2, and D in equation (11) are functions of some of the feedback variables 


in the core. 
The state vector of concentration fields of relevant fission products (like xenon, iodine and samarium) 


is represented here by x(t, r ) and its kinetics is represented by equation (3). 


The last three equations describe the dynamics of the feedback variables: thermal-hydraulics variables 
in the core (state vector w ), mechanical variables of the solids in the core (state vector y ) and other 


variables of interest outside the core (state vector w, ). 


Amongst the thermal hydraulics variables in the core we have the temperature field in fuel, moderator, 
reflector and coolant, as well as the flow velocity fields and void fractions of coolant in the core. 
Amongst the mechanical variables in the core we have the change of volume of the fuel pellets and 
the displacement, distortion and vibration of fuel pins, fuel assemblies, control rods, and in general 
structural elements in the core. 

Amongst the remaining variables of the system, for power and large production or research reactors we 
have those stemming from the primary circuit piping and circulating pumps and heat exchangers. 
Furthermore, in the case of a NPP, it could be necessary to add equations related with turbines and 
electric power generating machines with their electric loads, or a suitable subset of these equations, 
depending on the purpose. The details of each model depend strongly of the type of system and of the 
problem that is going to be considered. 

In any case, if the idea is to apply an analytical approach as far as possible, it is necessary to construct 
simplified mathematical models of the reactor core and the rest of the interrelated systems. So it is 
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advisable to lump parameters in as many state variables as possible to describe their evolution by ODEs 
instead of PDEs. The number of different nodes to be used depends of the frequencies of the transients 
that are going to be studied. (Akcasu et al., 1971; Stacey, 2018) 


Next we consider a steady state solution ¢, , Cro>%0> Yoo Wo» Wyo Of the equations, after removing the 
external neutron source S . This solution corresponds to a critical state of the nuclear reactor core. 
The problem to be considered now is the behavior of the transient solutions in some neighborhood of 


this steady state. This will be done in each of the three cases considered in this paper after making one 
or the other of the following approximations to nuclear reactor dynamics. 


2.2. Effective lifetime and constant production of delayed neutrons approximations 


Let us consider now two simplifications of the mathematical model of reactor dynamics. 


The average lifetime t, of delayed neutron emitters is of the order of ten seconds. If the time scale of 


the flux variations is at least an order of magnitude greater thant, , from equations (1) and (2) follows 


the effective lifetime approximation (ELA) done in equation (1) (Akcasu, Lellouche and Shotkin, 1971; 
Stacey, 2018): 


Gear 1H] & as.) = Wir] —L[ parson] 8 (12) 
u Ot Ot 


ELA will be employed in section 3 to study a static global bifurcation and in section 4 in an analysis of 
xenon oscillations. 


If the time scale of flux variations is at least an order of magnitude smaller thant, , the production of 


delayed neutrons can be considered constant (CDL). If ¢ is the flux in the steady state before a 


perturbation, equations (2) for the kinetics of delayed neutrons emitters become irrelevant in the CDL 
approximation and equation (1) for the neutron flux reduces to the following one: 


=F = (1p) M [dry ]-L[b.x, 0] + BM [dx Yom] (13) 


In equation (13) X,Y), Wy are the steady state values of the feedback fields in the core. 


CDL will be assumed in the study of the possibility of having fast power oscillations coupled with 
mechanical vibrations in the core of the reactor (the so-called mechanical kinetics effects in the above 
mentioned book by Marguet). 


2.3. Restricted nonlinear modal analysis 


The modal methods of solution of the field equations of nuclear reactor dynamics represent the fields 
as linear combinations of known space functions weighted by unknown time functions (mode 
amplitudes). This tentative solution is substituted in the evolution equations. Applying a suitable 
criterion, like the method of weighted residuals, a system of ordinary differential equations for the mode 
amplitudes is obtained (Stacey, 1967). Thus, a complex space- time field dynamics is reduced to the 
study of the evolution of a representative point in the space of mode amplitudes. 


In the cases that we are going to consider, the set of known space functions is a numerable and complete 
set of eigenfunctions of a linear operator. This operator, in general not necessarily self-adjoint, is 
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constructed from the field equations. The criterion that is applied to obtain the evolution equations for 
the modal amplitudes is the projection onto each eigenfunction of the adjoint operator (Henry, 1975). 


Often it is possible to work with a relatively small number of mode amplitudes, after reducing the 
original high order dynamics to a low order one. However, the linear terms almost always involve a 
certain degree of coupling between mode amplitudes. For example, if the Lambda — modes are chosen 
as the abovementioned set of space functions, then in the linear approximation to the evolution 
equations each mode amplitude that corresponds to the neutron flux is coupled with the others through 
a reactivity matrix (Henry, 1975; Suarez-Antola and Flores-Godoy, 2014(a); Stacey, 2018) 


However, there is a method of modal analysis, called restricted nonlinear modal analysis (RNMA) that 
can cope directly with the non-linear terms in the dynamics and often allows the determination of 
analytical formulae for the threshold amplitudes for instability associated with locally but not globally 
equilibrium points (Eckhaus, 1965; Denn, 1975; Suarez-Antola, 2005 (a)). 

For this reasons, a brief description of the RNMA method is introduced here, and subsequently applied 
to the examples analyzed in parts 3 and 4. 


To be able to apply this method, a proper choice of the linear operator and its eigenfunctions must be 
done so that the linear term in each one of the equations of evolution that corresponds to mode 
amplitudes, appears uncoupled from the other mode amplitudes. The importance of this uncoupling will 
be seen in the analysis of a reactor runaway developed in section 3. 

In the equation for the neutron flux field, the following linear operator is identified, fixing the feedback 
variables at their steady state values (Suarez-Antola, 2005(a)): 


A|¢]= al (1 — B)-M[9,x9.¥o.Wo|—u ; Lg, X91 ¥o>o | (14) 
Then, the reaction diffusion equation for the neutron flux can be written as follows: 


a 4 , K 
f= Ald]+ Nolo, x Xo ¥—Yoow w]tu YA, -c, tu-S (15) 


k=l 
Here NV, is a nonlinear operator. Fixing x —X),¥— Vy,W— W, , this operator is linear ing . 


It can be expanded as a sum of multi-linear operators of its arguments, beginning with a bilinear one. 
Usually a small number of multi-linear terms (two or three) will be enough. 


We pose the eigenfunction — eigenvalue problem for the operator A,. We use these eigenfunctions to 
develop the neutron flux ¢ and the fields of delayed neutron emitters c, ina series expansion in terms 


of the abovementioned eigenfunctions. 
Then we represent the fields x, y, w of feedback variables defined in the reactor core, also as a linear 


combination of suitable eigenfunctions with the corresponding time dependent mode amplitudes. 
Substituting this new ansatz in the evolution equations and projecting onto the eigenfunctions of the 
adjoint operator A, , we obtain a system of nonlinear ordinary differential equations. 


The regularity of the non-linear operator in the original field equation has as consequence that the 
differential equations for the mode amplitudes are of polynomial type (each second member is a 
polynomial in the mode amplitudes). 

In any case, we obtain a coupled system of ordinary differential equations with mode amplitudes of the 
neutron flux and the feedback fields x, y,w as state variables. In these equations a certain number of 


geometric, mechanical, thermal, hydraulic and neutronic parameters appear. Each combination of 
parameter values can be represented as a point in a parameters space. For each point in the parameters 
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space, there is a manifold of trajectories in the state space of mode amplitudes of the reactor. 
Depending of the problem, it is always possible to work with a finite and relatively small number of 
eigenfunctions, because in practice all physical systems behave as low pass filters, showing a significant 
damping of the higher modes (Stacey, 1967). This is what happens in the dynamic systems considered 
in sections 3 and 4 of this paper: after a fast transient from the initial conditions, only a few dominant 
modes (active modes) are enough to study stability and bifurcations. 


As a tule, the calculations needed to apply the asymptotic methods of restricted nonlinear modal 
analysis (RNMA) may be long and tedious. However, powerful symbol manipulation packages are now 
available to develop complex symbolic calculations in the computer. 

As these are asymptotic methods that depend of some sort of regularity in the non-linear operator of the 
field equations that allows a truncation of the system of coupled equations leaving polynomials in the 
mode amplitudes of low degrees. This regularity produces a kind of continuity of behavior between the 
linear and the non-linear regimes. The danger of this is that if the phenomena of interest are produced 
beyond the range of validity of the approximations, they could pass unnoticed (Denn, 1975). 

The other disadvantage is that if there aren't a few dominant modes the analytical approach probably 
will not succeed, and a numerical calculation should be undertaken. In this point it is necessary to assess 
if it is better to leave RNMA equations and use a computer code based on an ab-initio discretization of 
the field equations. 


2.4. Averaging method for nonlinear differential-integral equations with rectification in a 
neighborhood of a bifurcation point. 


Given a ROM with a small number of state variables, sometimes it is possible to make further 
simplifications, deriving a differential-delay equation in a single variable that summarizes the main 
characteristics of the dynamic behavior of the system. 

In the analysis of both, hard self-excited xenon oscillations and soft self-excited oscillations stemming 
from a possible coupling of nuclear, thermal and mechanical vibration processes, we will arrive to an 
equation of the following type: 


X(t) ~ F((t),x(¢)) (16) 


F(x,x)=—f(x)-%-@¢-g(x)-@;: | K(u)-g(t—u)-du (17) 


0 


Here f (x) is a damping function (not necessarily always positive), g(x) is a nonlinear restoring 


response with rectification properties (g(0) =0 but g(x) # —g(- x)) and K (¢) is a delay kernel that, 


perhaps after some oscillations, tends to zero when time tends to infinity. 
The frequency @, is always different from zero. In the Xenon model the frequency @, is zero but in 
the model of the mechanical kinetic effect it is different from zero. 


The damping function, the restoring function, the delay kernel and the frequencies are well defined 
functions of the parameters that appear in the mathematical model. 


However, for all possible combinations of parameter values, (x =0,x= 0) is always an equilibrium 
point that corresponds to a critical state of the nuclear reactor. 
When the dynamics given by equations (16) and (17) is such that x (t) oscillates with two widely 


separated time scales, the possibility of applying a method of averaging to study its asymptotic 
behavior must be considered. 
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Here we describe the method of averaging known as KBM (from Krilov, Bogoliuvov and Mitropolsky) 
to obtain approximate solutions in the form of nonlinear transient oscillations with a slowly (relative to 
the approximate period of oscillation) variable amplitude, combined with the energy method of Denn 
and Black to take into account the rectification properties of the nonlinear restoring response. The 
foundations of KBM method can be found in Minorsky (1983), Attlee-Jackson (1989), and Verhulst 
(1990). The method of Denn and Black can be found in Denn and Black (19733), Denn (1975). 


We begin with the following ansatz (tentative solution): 
x(t)=c(a(t))+a(t)-cosy(t) (18 a) 
y (t)= a -t+A(t) (18 b) 
We assume that both a (t) and 6(t) are slowly varying functions of time, their time scales being at 
least an order of magnitude less than = i 
0 
Furthermore, we introduce the restriction, characteristic of KBM method: 


2) _9,-a(0)-siny (0 (19) 


The offset term c(a ) in (18 a) is necessary to take due account of the rectification properties of the 
nonlinear restoring response g (x) . Often this is not mentioned in the expositions of the KBM method 


like the ones that can be found in the books of Minorsky and Attlee-Jackson cited above. 
We also assume that both, the duration of the interval of time [0,¢ Kl where the delay kernel K (t ) 


2 i 
differs significantly from zero and the period 7 = 2 are at least an order of magnitude less than the 
My 


order of magnitude of the interval of time during which the amplitude a(t) shows a significant variation. 
In this case, in the interval where the delay kernel is not negligible, a (t — t') & a(t), c(a(t = t')) & c(a (t)) 
, At—t') = A(t), so: 

yt-t)=y(t)—ay-t' (20a) 

x(t —1') = c(a(t))+ a(t)- cos{y(t)— a, -1'] (20b) 
In order to estimate c(a ) following the procedure suggested by Denn and Black, let us integrate both 
members of (16) in the interval (1,1 oe ) : 


[°" x(¢)-ar=a(¢+7)-2(¢)= [" Fe). a@)) ae 


From the equality f (x () 20 = “ i vd (x') dx'= La (x(t)), and interchanging the order of 


dt 
+T A 
integration in the integral of the delayed kernel, we find that the integral } Maar [x(t'), x(t')] -dt' can be 
t 


written as the sum of three terms: 


~(F(x(¢+7))-F(xl¢)), - 03” e(xle))-de and — 0 *K() [fall _ 1"). ar"|. dt 
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Now, making the approximations x(t +T ) = x(t ) and F’ (x(t +T )) =F (x(t )) , taking into account that 
) _ g(x(t—2"))- dt'= ) ‘ g(x(z'))-d¢' for anyt, using equation (18 a) we derive the following 
approximate relation: 


[~ e(c(a (t'))+ a(t')-cosy(t'))-dt'=0 (21) 


ae 
But a(t ) remains approximately constant in (1,1 zoho ) and dt'= —-dy so equation (21) can be recast 
0 


as follows: 
270 
li g(c(a)+a-cosy)-dy =0 (22) 
From this last restriction the function c = c(a) can be determined. 
When the amplitude a@ is small enough, approximating g(x) up to order three in a neighborhood of 


x=0, g(x)*x+g,-x’+g,-x°, we derive from equation (22): 
& §2 83 


Jy [(cla)+a-cosy)+2,-(c(a)+a-cosy) +25-(e(a)+a-cosy) |-dy=0 3) 


From (23) the following relation between the offset c and the amplitude of oscillation @ results: 


1 3 
tere tag 43-ge)+g-e = (24) 
When the amplitude of oscillation is small relative to 1, the offset is given by: 
1 
CET 8a (25) 
Once c = c(a) is known, it is possible to derive the following approximate equations for the amplitude 
aand the phasey : 
d 1] 27 A 
—az=- -| Flx,x|-siny-d (26a 
fos ae [Fle 4]-siny ay ) 
d QD 1 TS 
= ‘| Flx,x|-cosy-d 26b 
he ia ae [, Fle +]-cosy dy (26 b) 


In these equations x= c(a )+ a-cosy andx *—@,-a-siny 
Both the offset c and the amplitude of oscillation a are taken as constant for purposes of integration 
with respect to the phase of the nonlinear oscillator (Verhulst, 1990). 


From equation (26 b) we see that the local frequency “y = o(a ) is a function of the amplitude of 
t 


oscillation. 

The averaging method that allows for an offset that varies with the amplitude of oscillation seems to be 
not well known by most nuclear engineers and scientists. It was summarized here to ease its application 
in cases in which the nonlinear restoring term is not anti-symmetric. This last case is common in nuclear 
reactor dynamics, but not in nonlinear mechanics, where the analytical approach is often used. 

Here only the first approximation of a sequence of increasingly better ones was presented. As shown in 
the literature on averaging methods, successively better approximations can be obtained in the 
framework of KBM method (Minorsky, 1983; Attlee-Jackson, 1989; Verhulst, 1990). 
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However, the first approximation will be enough for our present purpose. The difficulties of the 
calculations increase fast with the order of the approximation, but as in RNMA case, powerful symbol 
manipulation packages are now available to develop complex symbolic calculations in the computer. 


Now, after having reviewed some of the mathematical tools necessary to achieve the objectives 
proposed for this work, we can begin with the first case to be considered: a static global saddle-node 
bifurcation that appears in the framework of a restricted nonlinear modal analysis of a reactor, 
subcritical at zero power and with a suitable nonlinear feedback. 


3. A static instability related with a global bifurcation 


Let us consider a bare reactor, whose extrapolated core fills a region B . Let us take into account the 
simplifications that allowed us to introduce equation (11) for an equivalent homogenized reactor. 


x f 
“and the 


4 
Taking into account the definitions of the effective multiplication factor k= 


she k=l desma 
corresponding reactivity p = ae , the prompt mean neutron generation time A = err and the 
U-V:Le 
f 


effective mean neutron generation time A, =A+/-t, = f#-t,, redefining a new field of external 


ie St,r 
neutron sources F’ (t, 7) = SWF) and introducing the parameter MZ —— , equation (11) 
A.-Vedy Keg 
can be simplified and re-written in the ELA general framework given by equation (12): 
O s s 4 “ 
A L67)=p-A07)+M?-V°OF)+A- FO?) Qn 


Now, let us suppose that there is a feedback mechanism whose main effect is a modification of the 
reactivity, so that its effects on the other parameters of the core can be neglected in a first approximation. 
Furthermore, we assume that the reactivity is a quadratic a polynomial function of the feedback variable 
Ww: 


P= Py ty W- Ay. W (28) 


Here p,, @,, and @,,, are positive parameters. As consequence, if w increases from zero, ¢ first 


wl 


increases and then decreases. If (, 2 0 the function p(w) takes positive values first, but when the 


feedback variable increases enough, the reactivity takes negative values. 


Let us assume that the field w(t Ae ) varies according to the following evolution equation, being 7,, a 


characteristic time scale of evolution of w, and the parameter K the link between the neutron flux and 
the rate of variation of w: 


La a ar ae 09) 


(aie 
Ot 


The boundary conditions for a bare reactor are O(t ,7,)= 0 and w(t 7,)= 0, for every instant f and 


every position vector 7, € OB (OB is the boundary of region B ). 


In the reactor model the feedback field w(t,7 ) could be, for example, a difference between a local 
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effective temperature and a reference temperature. In that case the parameter K could be given by 


= es) where c is an effective specific heat capacity, O is a coolant volumetric flow, u is 
. Uu . Cc . 


the neutron velocity in the one group description and « gives the heat power per neutron. The 
characteristic time could be given byT,, = e , being C an effective heat capacity of the core and h an 
effective heat transfer coefficient (Lewins, 1978). 


Now, our purpose will be to study the stability of the stationary solutions of equations (27) and (29). 
Let us suppose first that the source of external neutrons vanishes, so that ¢, (7) =0 and w, (7) =O isa 


possible steady-state solution (zero power solution). 

We introduce a suitable distributed source of external neutrons, during a while, to produce a perturbation 
in the field variables relative to the strictly zero-power solution. After the desired perturbation is 
produced, the source is removed. 


Let us consider what happens in the framework of nonlinear modal analysis. 


In this method the dominant mode amplitude is uncoupled from the others so that its threshold behavior 
can be studied in isolation. 


In this case the linear operator Ay [o| obtained fixing w at its steady state value w, (F ) =0 (for every 


position vector in the core) is: 


i M F 
A|6]=——V orig (30) 


The nonlinear operator N ; [o; w] now is given by: 
ie 1 1 2 
N,\9; w|=—@,,,-w:@-—-Q,,-w : 31 
oli wl= 7 aw O— 7a Wd (31) 


Then equation (27) may be cast as follows, in a form suitable to begin with a restricted nonlinear modal 
analysis: 


“ = Alé]+ Nlbsw]+F (32) 


In this case the eigenfunction-eigenvalue problem for the operator A, [¢]. with homogeneous boundary 


conditions, is equivalent to the eigenfunction-eigenvalue problem for the operator — Vv’ [o| with the 
same boundary conditions for the neutron flux. 


The eigenfunctions ¢, (7 ) of —V” [o| verify the Helmholtz equation— V9, (F ) =U, (F ) 
The positive eigenvalues Le can be ordered in an unbounded increasing sequence 0 < uy) < LG <.... 


() The eigenfunctions are orthogonal and can be normalized: (Q,+Pn) = } Q, (7 ): Q,, (F )- dV =6,,, 
B 


being 6,,, Kronecker’s 0. 


These eigenfunctions are the same for both operators — Vv’ [d| and Ay [?] , 
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; —M? uP 
To each eigenvalue ue of the operator — Ve [o| corresponds an eigenvalue ame of the 


e 
operator A, [0] As a consequence, to study the stability of the zero-power solution, we represent the 


fields as series of eigenfunctions @, (7) with time varying amplitudes: 


o(t.7)= > 4,(0-0.7) (33 a) 
W(.7)=¥B,(0-9,(F) (33 b) 


n=l 


Substituting the ansatz (33) in equations (27) and (29), projecting onto each eigenfunction @,,, the 


following infinite system of nonlinear ordinary differential equation is obtained, for p =1,2,3... 


dA “ a 
A : 2 =-(M?y? - p,)-A, +a: Yo Limn An? By = Ey" > Ligne Bs B, + F,(t) (34a) 


dt m,n=0 m,n,q=0 
OB, 
ae es =K-A,-B, 
(34b) 
Lan = [Pp “Pn Pq AV (34c) 
B 
Lnng = | Pp “Pn * Pn Py AV (34d) 
B 
F,= [¢, ()-F(t,r)-aV (34e) 
B 
The integrals 7, and /,,,,,,, , temain invariant under any permutation of the indexes. 


If -A, (0) =0 for every mode index p =1,2,3,....and at least one of the forcing terms F , #0, then a 
perturbation will be produced from the zero-power solution. 

Once every projection F, (t) has vanished, the perturbation has already been produced and the next task 
is to determine its future evolution. 

To continue, let us suppose that at zero-power the reactor is sub-critical. This means that the 
inequality p, <M : We is verified, and because Le, > rif if p #1, all the coefficients of the linear 
terms are negative in Equation (34 a). 


Furthermore, M7? LU — Py increases monotonically and without bound with the index p . 


Now we suppose that the outer time scale of the neutron flux 7, = A, if (uw Lh po) is an order of 


magnitude greater than the time constant of evolution of the feedback variable. 

Applying the second principle to link time scales enunciated in subsection 1.2 of the introduction, we 
can consider that w is in equilibrium with the neutron flux. 

With this assumption, the system of equations (34) reduces to the following: 
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A 00 0 
A, ee (ui Solid: +Q, “KG yee 1A “A, ~ Ay» -K? 3 nies ‘Ay “Al, ‘A, + F,(t) 


m,n=l myn,g=l 

(35) 
Let us take as our time origin the instant in which all the external forcing terms F’, (t) vanish. 
Then, as the linear terms have negative coefficients, for a perturbation near enough to zero the state of 
the system given by the mode amplitudes 4, (t) will return to the origin of the space of mode 


amplitudes. 


Now, we suppose that Mf” - L :— Pp 1S very near zero, so that in a neighborhood of the equilibrium the 


amplitude of the first mode shows a critical slowing down. 
However, the coefficients of the linear parts of the equations for the other mode amplitudes are not 


2 2 3 3 : 
necessarily near zero, because M° - 1 » increase with p and the order of magnitude of M aa Lo- Po 


will be greater than the order of magnitude of M? - °— fp if this last one is small enough. So, we can 


assume that there is only one mode amplitude that shows a critical slowing down. 


Let us suppose also that the external source excites mainly the zero-order mode amplitude A (¢ ) ; 


because F(t) is much greater than F, (t) for p #1. Then we can expect that the amplitude A (t) will 


be dominant, and it will be possible to uncouple it from the others mode amplitudes. (Eckhaus, 1965; 
Denn, 1975; Haken, 2003) 
We thus obtain the evolution equation for the uncoupled dominant mode: 


d 
9) at = (Mr? ~py)-A, Og TA Oi KR ye A, + F(t) (36) 
(a > OF > 0) 
We can expect that the others mode amplitudes evolve slaved by the dominant mode: 


d 
A, <4, = (M?u? —py)-A, + Oy) °K Ty, ° A? Oy K? 1, A +E) p¥l (7) 
The others nonlinear terms in these equations are negligible (Eckhaus, 1965; Denn, 1975; Haken, 
2003). 


So, if the external source was active fort < 0, and is withdrawn in the instant t = 0, we have the initial 
conditions A, (0), A, (0) ,--. With | A, (0)| much greater than the other| A, (0) ; 


Once the external source disappears from the core, the uncoupled dominant mode amplitude evolves 
according to the homogeneous equation, which is the normal form of a static global bifurcation: 


d 
iN, a = (M713 ~py)-A Bag RT e A, Sage oD (38) 
The right-hand member of the equality has at least one real root A,=0. It is the only real root if 
ar 
M? + — p, >“ . It is always stable. 
4 yy Li 


If it is the only real root, then lim A, (t) =( and the higher mode amplitudes will be forced to approach 
t+too 


to zero, as is easy to see from their evolution equations. 
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2 72 
But if M? +? —p,<—“ there are two other real roots, A,, and A, ,. 
4: A.D 
Ay, is unstable and A, s 18 stable. The unstable mode amplitude is given by the formula 
Ae Ali) — ade 40,0, (Meu? ~~) (39) 
lu 
2K QL) 


As we supposed that M7? - 41;— p, is near zero enough so that the order of magnitude of M *, | Po 
with p #1 will be greater than the order of magnitude of M’ - 4;— pp. 

As consequence, the time scales 7, = A, / (u ? Lt, - Po) of the higher order mode amplitudes (p #1) 
are an order of magnitude smaller than the time scale 7, of the dominant mode amplitude A (t ) . Then, 


after a short transient relative to 7, , the slaved mode amplitudes will behave approximately as: 
A jer AP (t)= °K? Lyi A(t) 
‘ (u < H, 7 Po) 


pill ; 


(40) 


Equations (36) and (40) determine a slow manifold where the relevant part of the dynamic of the 
original system, given by equations (35), can be analyzed. 


If A, (0) <A_,,,, then A, (t) will return to zero dragging the other mode amplitudes to zero. 


But if A, (0) > Ay» A (t ) will grow, approaching to A,, and dragging the other mode amplitudes to 


follow it: the whole neutron flux will leave the set of states attracted by the zero flux condition, and will 
tend to another stable steady-state flux. 


If we take yu; as a parameter to vary, it is inversely proportional to the square of characteristic linear 
dimension of the reactor core. If the core is small enough, yz; will be big enough so that 


2 
"ee be 


wl 111 


M? - i - py > . Then the zero flux will be the only steady state solution, and it will 
1 a eae 
w2 1111 


attract all the other states. So, no matter the amplitude of a perturbation, it will never cause instability. 


But if the core is big enough to reverse the above inequality, a threshold amplitude appears and with it, 
an instability to finite perturbations is produced, even if the zero flux remains locally stable. 
Thus, the restricted nonlinear modal analysis method allows us to find and quantify a static global 


bifurcation (a saddle-node bifurcation) with two branches, one corresponding to A,,, and the other to 


Aj 3 


Lg - 
When the bifurcation parameter P = — is lower than a critical value P. = 


a) > 
Hi (> nm tiny 
0 
4 Oyo Lin 


is a single equilibrium point: the nuclear reactor off, corresponding to A,=0, and all the other mode 


amplitudes zero also. 
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But when the bifurcation parameter attains its critical value, an unstable (threshold) equilibrium modal 
amplitude and a stable equilibrium modal amplitude bifurcate and behave as shown in Figure 1 as the 
bifurcation parameter P grows. 


5.00 =] 
1 Al 
4.00 — 
3.00 
2.00 
1.00 — 
' 
' 
1 A = 0 
1 ee 1 
ar, 
0.00 2.00 4.00 6.00 8.00 P 10.00 


Fig. 1: Bifurcation diagram for equilibrium amplitudes. The unstable branch of the saddle-node 
bifurcation is given by A , as function of P. The stable branch is given by A , as function of P. 


The bifurcation diagram of Figure 1 corresponds to a slab reactor of thickness /. For this reactor the 


2 2 
; 1 - ; ; ; l 
eigenvalues are #6 = = so the bifurcation parameter is P=—. 
a 
a2) 


The eigenfunctions for the slab reactor are cosines or sines of = -x for p odd or even, respectively 
(Duderstadt and Hamilton, 1976). In this case x represents the distance measured form the middle of 


the slab, so that oe <x< : : 
2 2 


Now, an example of digital stimulation of interactive mode dynamics will be considered to check the 
assumed dominance of the first mode amplitude in a nonlinear systems of interactive mode amplitudes. 
Let us observe first that the following version of the neutron balance equation in space-time is derived 
for a slab reactor, from equations (27) (28) and (29), with the feedback variable in equilibrium with the 


flux, and with both members divided by | Po| : 


A, ap(t.x) Ba ats Se) ge (o)-{ Sa) 4 


|Po) a(t) 7 [Po Po| 


(41) 
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When the parameter (, takes negative values, Pit =-—1 and equation (41) is formally identical to an 


‘|Po] 
equation that describes the threshold effects in a mathematical model of excitation of nerve fibers by 
external electrodes. (Suarez-Antola and Sicardi-Schifino, 1996) 
Mode decomposition was done for this equation, and the interactive mode amplitudes were truncated 
after p =5. The details of the simulation results are summarized elsewhere. (Suarez-Antola, 1997) 


For the spatial distribution of the forcing term the following formula was employed: 
F (teos{ = “J+A ( Jeos{ 2 ere Jeos{=2=) 


The system starts from rest (reactor turned off, membrane at its rest state in the nerve case) and the 


functions F, (t), EF, (¢) and F; (t) are chosen as step functions equal to the constant values F;°, F;’ 


and FY respectively during an interval of time of duration — , after which the forcing term (external 


ial 
neutron sources in the reactor case, external imposed electric current in the nerve case) is withdrawn 
and the systems evolves unforced. 

Due to the symmetry of the distribution of the forcing term, only the first three mode amplitudes of odd 


index (A,, A, and A,) are dragged from their zero value. Both A, and A, remain undisturbed. 


A negative value of the parameter (, strengthens the stability of the rest state (zero-power equilibrium 


Po 
in the reactor case) and is necessary in the nerve excitation case. A consequence of having ——~ = —1 


[Po] 
is that there is not critical slowing down of the first mode amplitude unlike what happens when the 
parameter (, takes positive values and M -s bie Pp 1S very near zero. 


However, the simulation results show, as will be seen below, that even in this case presumably 
unfavorable, the dominance of the first mode persists. 

°K. @yyK° L gle L? 
The other parameters were selected such that —4 =], = =l, 


|Po| [Po Taal 
Figure 2 shows the results of the simulation for F° < FY < FY, being F° an order of magnitude less 
than F.’. 


= 


mode amplitude 
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Fig. 2 An example of temporal evolution of mode amplitudes. Dimensionless time appears on the 
abscissa and mode amplitudes in ordinates. 


The left part of the figure shows the evolution of the coupled mode amplitudes when the system is being 
forced by the external source of neutrons. The first mode amplitude attains a value smaller than the fifth 
one, and this attains a value smaller than the third one. 

The right of the figure shows the evolution of the coupled mode amplitudes when the system is free 
from the external source of neutrons. It can be seen in the figure that the first mode amplitude growths 
towards a new equilibrium value and behaves as dominant over the other two, that tend to their own 
equilibrium values. 

In the following section we apply the method of restricted nonlinear modal analysis to analyze xenon 
oscillations. 


4. Hard self-excitation of Xenon oscillations revisited 


The effects of Xenon are a potential source of instability in thermal reactors. Xenon oscillations are 
produced by the delay between xenon burn-up and xenon build up from iodine decay when a change in 
local neutron flux produces an imbalance between both processes. This can establish an oscillatory 
regime in reactor power with dangerous peak values, mainly when the reactor is used for load following, 
with frequent power changes. 

Two types of in phase self-excited Xenon oscillations have been experimentally found, described and 
analyzed by means mathematical models of point kinetics with feedback: soft and hard oscillations 
(Chernick, et al, 1961; Rizwan-uddin, 1995). 

In the framework of point kinetics with feedback, soft oscillations appear when a change in the 
parameters of the reactor drags a globally stable steady state to the boundary of stability in the plane of 
parameters determined by , (the zero flux and zero Xenon poisoning reactivity) and y (a prompt 


power reactivity coefficient). After crossing the threshold barrier in a suitable place, the steady state 
becomes unstable and a stable limit cycle of oscillation appears whose amplitude grows from zero: a 
supercritical Hopf bifurcation is produced. 

Hard oscillations appear when an unstable steady state crosses the threshold boundary at another place 
of the p,-y plane, becomes locally stable and at the same time an unstable limit cycle is born, whose 


amplitude also increases from zero, but disappears when the steady state is inside the local stability 
region far enough from the stability boundary: a subcritical Hopf bifurcation is produced. 

Besides, several types of out of phase self-excited Xenon oscillations (in axial, radial and azimuthal 
directions) have been experimentally found, described and analyzed using space-time mathematical 
models and modal analysis (Bell and Glasstone, 1970; Lewins, 1978) or multipoint kinetic equations 
(Chakraborty et al., 2018). 

In this section our focus will be in the study of only some aspects of Xenon oscillations. 

In 4.1, we consider the instability thresholds to finite amplitude perturbations for in phase Xenon 
oscillations by the methods of subsections 2.3 and 2.4. Only a brief consideration will be given to the 
thresholds and periods of out of phase oscillations in 4.2. 


To study stability and bifurcations related with Xenon oscillations by analytical methods, let us consider 
again a bare reactor, whose extrapolated core fills a region B . As already said in section 3, the boundary 
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conditions for a bare reactor are g(t ie) =0 and w(t 7, )= 0, for every instant ¢ and every position 
vector 7, € OB (OB is the boundary of region B ). 


A possible set of equations suitable to study Xenon oscillations in an equivalent homogenized reactor 
are, making the ELA approximation that appears in subsection 2.2 (Suarez-Antola, 2005(b)): 


A, (F)= p07) HM? -V°A(6F) (420) 
Oi 
a Sty Ly p-A,-i Ge) 
Sag. Dt A inAex-o 0-9 (42c) 
OAT AT. 
pe Rg AT. 2d 
es F d (42d) 


In the balance equations for Iodine and Xenon, x is the concentration of Xenon, and i is the 
concentration of Iodine. @, and @, are the Iodine and direct Xenon yields and 2, ,A, are the 


corresponding decay constants. In the equation that gives the evolution of the difference AT between 
the local effective temperature and a reference temperature, 7 is temperature’s time constant (that 
represents the time lag in the temperature feedback produced by the thermal capacity of the reactor). 
AT, and ¢, are an equilibrium constant temperature and neutron flux respectively, while o, is 


x 
Xenon’s microscopic absorption cross-section. and @ is a temperature feedback coefficient. 
The local reactivity is a function of Xenon concentration and of temperature: 


o.-x a-AT 
AT |= x 
p ee | Po * = F S. 


The meaning of M* and v> ; is the same as in equation (27) of section 3, and —a@ is a feedback 


(42e) 


temperature coefficient (with @ positive to have a negative feedback). 

The reactivity when x =0 and AT =Oisp, >0. 

The problem now is to study the stability of the steady-state solutions of these dynamic equations. Given 
the homogeneous boundary condition applied to the neutron flux, the zero solution is always an 
equilibrium solution. But if 9, is positive and the dimensions of the core are big enough, there is always 


a unique positive solution g, (r) i (r), Li (r) and AT, (r) that represents the corresponding just- 


critical state of the reactor. The neutron evolution equation is re-written as follows: 
ae ‘ 
A, 2 = Ay [d]+ No[bix—x,AT-AT,] (43a) 
t 
Equation (43a) is a particular case of equation (15) of subsection 2.3 (Restricted nonlinear modal 
analysis). Here the linear operator (self-adjoint) and the nonlinear one are, respectively: 


A,|p]=M5V°6 + p[x,,AT,]-¢ (43b) 
A” oO a 
= AP SAT. 
Ny [osx Xp AT AT | vy, (x xy) p 7, | 0) p (43c) 


The idea now is to consider the complete set of eigenfunctions ‘Y, (and their eigenvalues @,, ) of the 


three-dimensional Sturm-Liouville problem A, [w]= o-Y with Y defined in B and with zero 
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boundary conditions, and represent the neutron flux, the concentrations of Xenon and Iodine, and the 
local temperature, as follows: 


b(t.r)=&(r)+d.®, (t)¥, (7) (Ada) 


x(t,r) =Xy (r)+ ya (t)Y,, (r) (44b) 
i(t,r) =i,(r)+ 3 i, (0)¥, (r) (A4c) 
AT(,r)= AT, (r)+ >a (1) ¥, (r) (44d) 


The eigenfunctions are orthogonal and can be normalized: [¥, (r) oe (r). dV=6,, 
B 


Substituting the ansatz (44) in the dynamic equations (42) and projecting onto each eigenfunction Y’, (r) 


, an infinite set of nonlinear ordinary differential equations for the time dependent mode amplitudes is 
obtained. 


The operator Ay [¢] was chosen such that the mode amplitudes D (t) for the neutron flux appear 


uncoupled to the first order (that is, uncoupled to linear terms). 
In the modal expansion, the eigenfunctions are ordered according to the magnitude of their respective 
eigenvalues. Then, the expansion is truncated at a certain eigenfunction Y,, (r), using the following 


criterion. The higher harmonics will be excited only by reactivity perturbations localized in a region 
small enough, so that the reciprocal of a representative dimension of this region is longer than the 
eigenvalue separation of these harmonics (Zauderer, 2006). As consequence, modal expansions with 
only a few eigenfunctions may be almost as accurate as the output of full three-dimensional modal 
methods, if the initial reactivity perturbation is not too localized. For this reason, low order modal 
expansions may be enough to study common xenon instability problems, but are no adequate to describe 
a transient produced by, for example, a sudden drop of a control bar. 

After truncation above the order M, we obtain a finite system of nonlinear ordinary differential 
equations to describe the evolution of the mode amplitudes in a finite dimensional Euclidean space. 
From now on, we take the zero solution (zero power state) as reference. 

Of course, this is not the best thing to do in order to study instabilities during the operation of the reactor. 
However, it has two advantages. First of all, we obtain a set of equations that, in the case of gross xenon 
oscillations can be easily compared with the equations already studied by others. Second, we can see 
how it is possible to obtain other steady state solutions, different from the origin, and study their stability 
using the same modal framework for all of them. 


To fix ideas, let us consider the case of a PWR, withA, ~0.ls, andz = 3s. This time scale of 


temperature is much smaller than the time scale of iodine and the local time scale of xenon (Henry, 
1975) 

However, the local time scale of the neutron flux, despite the value of the effective generation time, 
may be much larger thant . 


2 72 
It is of the order of the time s/o + plx, ar} : 
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If the state of the reactor is near an attainable steady-state, the denominator may be fairly small because 
both the buckling and the effective reactivity have opposite signs and almost the same absolute values. 
Nevertheless, if the state of the reactor is far enough from a steady-state, the time scale of evolution of 
the flux may be fairly small. In any case, we are going to suppose that the local temperature is relaxed 


* 


to its equilibrium value AT,, = ¢. Other authors have done this approximation (Chernick et al., 


1961; Rizwan-uddin, 1995) 


* 


J a 
Introducing the new parameter y = 
V 
f T* 


and applying the procedure described in subsection 2.3, 
and taking the first two eigenfunctions only, we obtain for the mode amplitudes: 


For the fundamental (or global) mode: 


d® ae oe O*°. 
oe {6 —-S lh 8,=78ay ®, } )— F887 03 (45a) 


“dt Vidas vd, 
di : 
ae = @ De ‘Dy — A; +i (45b) 
aX, . 
ae =P, Dar Dy +A, Ig — A,X — F,Pyog * Xq* Py — Fy °° P; (45c) 


Notice the coupling between the fundamental and the first harmonic in the equations of evolution for 
the modal amplitudes corresponding to neutron flux and xenon concentration. 


For the first harmonic (or regional mode): 


d® ¢! Oe om 
A eS | Pi Og i S20 Bg WG | Deir HD (46a) 
dt vy voy 
ey eee ee er (46b) 
Ot 
= = Q, ya D, +44, — A,X, — 0,01 (x,®, a xX), ) (46c) 
By definition: 
A xmn = ie (r)-,,(r)-,(7)-dV for p,m,n€ {0,1} (47) 
B 


From its definition we see that 9 remains invariant under any permutation of its sub-indexes. 


pmn 


Observe the nonlinear coupling between the fundamental and the first harmonic in the equations of 
evolution for the mode amplitudes corresponding to neutron flux and xenon concentration. In this 


equations p% = ~)—-M>-u5 and p? = p,—-M>-u), where 5 and yy are the fundamental and 
the first harmonic eigenvalues of —V* = uu" Y with Y =0 on the boundary. 


The relation between one of this eigenvalues, 7? and the corresponding eigenvalue @,, of the operator 
A, introduced in equation (43b) above, is @, =? = )—-M>-u-. 


p,. is the so called “static reactivity” of mode n, without feedback. 
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‘YY, is symmetric and positive and ‘Y, is anti-symmetric considering a mid-plane through the core, so 
that Ayu) >9, Ao, =0, A, > 0, A, =9. 

Using only two eigenfunctions to approximate the fields in the reactor’s core, we represent the neutron 
flux by QW, (t )P, (r)+ ®,(t ee, (r) and iodine’s and xenon’s_ concentrations by 
ii (t ). a (r) ae (t ). YT; (r) and X, (t ). ey (r) Py (t ). we, (r) respectively. At this level of description, in- 
phase xenon oscillations appear mainly in relation with the amplitudes of the fundamental or global 
mode O, (t ), Rae (t ), ke (t ) , while out-of-phase oscillations appear mainly related with the amplitudes of 
the first harmonic ®, (t ), x, (t ), hi (t ) or regional mode. 


As a first approximation we may assume that in phase oscillations only the global mode is excited, and 
in out-of-phase oscillations the global mode remains near its steady state and only the regional mode 
(the first harmonic) is excited. 


4.1 Global Xenon oscillations revisited: an instability threshold to finite amplitude perturbations 
If we neglect the nonlinear coupling terms in the global mode equations (45) we obtain a simplified set 
of equations. Instead of (45a), (45b) and (45c) we have: 


d® ee. 1 
A, : a= \ Foo "Xo = FO ‘D, -D, (48a) 
dt vo, 
di, ; 
dx, 
Lap. = @Q, Dig ®, a Aig = A,Xq ra F,AyooXpPo (48c) 


These modal equations are mathematically homologous to the point kinetic equations posed by 
Chernick et al. (1961) to describe in phase xenon oscillations and studied more completely latter by 
Rizwan-uddin (1995). During their above-mentioned work, Chernick et al. made a drastic 
simplification: they made A, =0 and derived a second order nonlinear oscillator differential equation 
in terms of Xenon concentration. They didn’t solve it, neither with approximate analytical methods nor 
numerically. However, they studied the curve of threshold stability in the plane of parameters , and 
y predicted from zeroing the dissipation term in the oscillator equation. They compared this prediction 


with the curve of threshold stability derived from the numerical analysis of the full set of point kinetic 
equations. They found that the predicted curve was compatible with the curve obtained for the full set 


of point kinetic equations when the reactivity feedback coefficient was large enough (circa 10°) and 
the prompt power feedback coefficient was small enough (approximately 0.2x10°'°cm’ xs , see also 
Rizwan-uddin,1995). 

In equations (48) the homologous parameters to the point kinetic p,and y are, respectively, 
2 = P) -M;- us, and 7. So, for large enough values of ) and small enough y it would be expected 


that the same approximation made by Chernick et al. could be made also in the case of equations (48). 
From (48a) and A, =0 we derive an approximate equality to express the global mode amplitude of the 


flux as a function of the global mode amplitude of the xenon concentration: 


1 -.. 30% 
D, = "| Po - ooo "Xo (49) 
Y Ao voy 
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This approximation can be understood from the perspective that provides the theory of the slow 
manifolds summarized in 1.2.3, if D, behaves as the fast variable z and the concentrations of Xenon 


and Iodine behave as the slow variable y. 


Then G — FF 1X6 = VO seo o,}, =0 correspond to F(z,y;c)=0 and equation (49) 
f 


O 
correspond to one of the branches z = h( yc) To the local stability operator aot (h ( y3c) ; y;c) it 
Z 


corresponds the function —7 - Ay) -P ke so, if 7 is positive, the branch given by equation (49) is a stable 


slow manifold in the sense that was explained in 1.2.3. 


Substituting (49) in (48b) and (48c) and eliminating the global mode amplitude of iodine from the 
resulting two equations the following equation for the global mode amplitude of xenon concentration 
is derived: 


ae +(a,=6.-%9) S24 03 (x9 —b, xd.) =0 (50) 


This equation for nonlinear Xenon oscillations is based on a modal analysis of a spatiotemporal model 
of the reactor core, so parameters that appear in this case differ from parameters that appear in the 
nonlinear Xenon oscillation equation derived by Chernick et al. (1961). 


All the parameters in equation (50) are positive. With p= @,+@, being the sum of Iodine and direct 


Xenon yields, the parameters are given by the following equalities: 


Y V 
eck ae 
v> 000 
b, = F Z (51b) 
[A-r+(2+08)] 
om 4 
2 
pie Op. (Sle) 
Vy 
oe = 8.4, [4 r+(2+08) (51d) 
om Vv 
be }e 
Fee (51e) 
«-y+(2+p6) 
om Vv 
A,-o-X 
Then: ood, = [Aes (52a) 
Yq 
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[2-0 
ag 2 (52b) 


The equilibrium Xenon concentration X, of the nonlinear system (50) must be a root of the equation 
1 
X, —b, -x¢-d, =0. This equation has the roots X= aE +./1—4-b,-d, 


From (49) it follows that the static reactivity 9) of the global mode is the sum of two terms, always 


positive: (5 = ———- Ayo *X + YA *Po So p is positive. If this static reactivity decreases tending 


V Day 
to zero, both x, and ®, must approach zero. According to equation (51le) the parameter d_, tends to zero 
with the static reactivity, but according with equation (51b) the parameter b, decreases taking positive 
values when ; increases. 


Then, the only root that tends to zero when 9} tends to zero is: 


H=5 5 (vi #b-4.) (53) 


This is the equilibrium value of the global mode of xenon. From equation (52b) and choosing the values 


A 
of the parameters that appear in Chernick et al (1961) to estimate “e and —--y we find that b,-d, is 
Vv 


Oo 


x 


at least one order of magnitude less thanl. Then the value of X, must be always near the value of the 
parameterd, and b,-X, must be small relative to 1. 
If we put x, (t) =X, (1 +& (t)) in equation (50) we obtain: 
dé 3 = dé = = 
ae +((a, -¢, -%)—-(¢, Fy) €) + wo ((1-2-b, ‘¥))E-(b,-%)-E7)=0 (54) 
The linear approximation in a neighborhood of €=0 is a linear oscillator of natural frequency 
@ = @,9'1—-2-b,-xX, and damping coefficient a, —C, +X: 
We NE of = 
—~ +(a, —C,°X))-— + @5:(1-2:-5,-x,)-F =0 
dt’ ( 0) dt 0 ( 0) S 
If the static reactivity of the global mode is small enough (107 or less), c, -%) =c,-d, is less thana, 


and the damping coefficient must be positive. 
Let us suppose that the damping coefficient is positive, so that the steady state of the reactor, given in 


dé(t) 
dt 


Now we introduce the following functions and parameters, to be able to use the averaging approach 


this case by the variables € (¢ ) and , is stable to infinitesimal perturbations. 


reviewed in subsection 2.4: 
f (E)=(4-4-%)-(e -%)-E (55a) 
(€)=(€-2,-€7) (55b) 
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= (2. %y) 
7 55 
oP a0) a 
, ace des 
Then equation (54) can be written this way: a + (Gc) a +052 (é) =) 
So, it can be recast in the framework of equations (16) and (17) of subsection 2.4: 
d° a( d d os 
a? (S25 f (6) 6-@y 8(6) (56) 


Given the ansatz & (t) =c+a-cosy, from (22) we have a restriction to find the offset as a function 


of the amplitude of oscillation. From (24) and (55b) we derive: 
1 
Ct Rpt 58s a = 0 (57a) 
If the amplitude is small relative to 1: 


-Z,-4 (57b) 


2 
bs SO Z 
From (26 a), (26b), (56) and the restriction a is (t) =—@):a-siny we obtain: 
t 


d _ 1 a f (c(a)+a-cosy)-(-@,-a-siny) | 3 58 
dt 2n-d * +0 g(c(a)+a-cosy) — i 
ae 1 f(c(a)+a-cosy)-(-@,-a-siny) | F Ga 
Baan sr “COSY 

dt 2 2n-@-a  |4+0%-g(c(a)+a-cosy) - 


From equations (58) and (55) it follows: 


Bt Ute ave oa jy| Se (ence) (ote See 
qos 5 [ (4. Cc, xe) (2 cae c(a) |= 5 (a C, x; (e; cae 5 ga (59a) 


d s fe - 
o(a)=T y= @ (1-2-8, e(a)) =a, (I-%-a’) (59b) 
Equation (59a) can be recast as follows (the normal form for a subcritical Hopf bifurcation): 


a= -K-a-[a?—a"] (60a) 


In (60) the new parameters «(a new kinetic parameter) and a, (a non-zero equilibrium value of the 


variable € representative of global mode amplitude) are given by the formulae: 


Pa (¢.-%0)- 8 (60b) 


(60c) 
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Let us remember now that we defined the dimensionless variable & (t) through the formula 
Xo (t) = 05 (1 +& (¢ )) and introduced its offset c , amplitude of oscillation a and phase w through 
the formula (1) =c+a-cosy. As consequence, x, (t) =X ° [1 + c(a (t)) +a(t)-cosy (¢) | 


The modal amplitude of the global mode of the flux is obtained substituting this formula for x, (t ) in 
equation (49): 


®,(t)= u G 4 [-el a(t) ra()<060(0)] (61) 


Y¥ Poo ve f 
The approximate behavior of the global mode of the neutron flux (and thus the approximate behavior 
of global mode of the reactor power) are determined by equations (57b) and (60). 
From (60a) we find that if an initial value of the amplitude a (0) is less thana,, then the amplitude of 


oscillation decreases monotonically towards zero, and its frequency @(a) increases monotonically 
tending towards @, . 


If a(0) =a, then the amplitude of oscillation remains constant a(t ) =a, as well as its frequency 
o(a,) = & (1-8 -a.). 
If the initial value a(0) greater thana,., the amplitude of oscillation increases, and its frequency 


decreases monotonically as time passes. 

Thus, the average method predicts the existence of an unstable limit cycle in the global mode amplitudes 
of Xenon and neutron flux (and as consequence, of reactor power). 

If a perturbation from the steady state leaves the state of the reactor inside the circle, it will return to the 
origin. But if the perturbation leaves the state of the reactor in the region outside the circle, the state 
will move farther away: the uncoupled global mode of xenon oscillation and neutron flux, under the 
established conditions, presents a subcritical Hopf bifurcation. 


Furthermore, we have derived closed analytical formulae for the threshold amplitude a. , the kinetic 


parameter and the variable frequency of oscillation @(a) that can be used to predict the time 


behavior of the global mode amplitudes of xenon, neutron flux and reactor power. 

When the amplitude of oscillation grows sufficiently, the simplifying assumptions that were made to 
be able to deduce the analytical formulas (57), (60) and (61) for global xenon oscillations, cease to be 
valid approximations and it is necessary to resort to numerical methods. A Bautin bifurcation was 
discovered by digital simulation. The bifurcation software and the digital simulation of the dynamics 
shows that a stable limit cycle surrounds the unstable limit cycle and the radii of both approach each 
other until they become equal. After this cycle coalescence the equilibrium point of the reactor becomes 
globally asymptotically stable (Rizwan-uddin, 1995). 


4.2 Out-of-phase xenon oscillations: instability threshold and period. 
To attain a certain degree of completeness in the proposed modal analysis, let us consider briefly the 


excitation of the regional mode while the modal amplitudes of the global mode ®, and x, remain at 


their steady-state values D, and x,. Then, equations (46) become: 
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d® e. Oe = = om = 
A, =| Pf -—- + y91 -%q— 2+ Y -Fo1 “Do |- Dy —-—— Fo, "Dy (62a) 
dt vy voy 
di 
- = @; par @,-A, +i, (62b) 
dx, = 8 
OE = 9, Lip, +4, +5, - A,X, -F,4 (x *D, +X -®,) (62c) 


As we did in the case of global xenon oscillations, we will assume that the regional modal amplitude of 
the neutron flux is relaxed to equilibrium with the regional modal amplitude of xenon concentration: 


(=| 5.) [(au[ $5 +7-5,)-oF}} nt (63) 
f 000 


If we substitute (63) in equations (62) we obtain a system of linear differential equations in the 
amplitudes corresponding to the regional mode of xenon and iodine concentrations. 
Now, following a procedure like the one employed in 4.1 and making a linear approximation we find 
the equation for a damped harmonic oscillator: 
d ca 
dt 2 
The linear damping coefficient 2-¢-@, is given by 


2-6-@°: Sia: x, =0 (64) 


yoo 


This coefficient may change its sign if the temperature feedback y is weak enough and the flux is high 


220, =, +4, +0.0y 8] (209,75 % 01} /O9( 247-8, | (65c) 


enough, with y-®, small enough. 


The combinations of values of the parameters of the reactor such that the linear damping coefficient is 
zero defines a stability boundary for infinitesimal perturbations. 
The natural frequency is given by: 


O; =A, (A, +0.) (2278 -®, —p; Gs [Bia }-o ] (66) 
000 


Here p=, +@, is the sum of the iodine yield and the direct xenon yield. 


The true frequency of the oscillator and its period may be calculated by well-known formulae for a 
damped linear oscillator. 

A more complete study of out-of-phase xenon oscillations must be done by bifurcation codes 
(Kuznetsov,1998) and digital simulation of the dynamics in the framework of a suitable ROM 
(Chakraborty et al, 2018). 


5. Mechanical kinetic effects and soft self-excitations of power oscillations 


In his recent book of nuclear reactor physics, Marguet poses the problem of investigating the possible 
dynamic effects of a variation in geometry due to the variable mechanical stresses in the solids present 
in the reactor core (Marguet, 2017). This is a different kind of problem that the well-studied one that 
arises in connection with the variations in the void fraction in the BWR coolant. 

The slow and stable power oscillations (period circa 30 s) observed in the EBR-I when steady state 
power was high enough, apparently was produced by the combination of a prompt positive fuel 
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feedback coefficient (due to the bowing of the fuel rods toward the center of the reactor) with a larger 
delayed negative feedback (due to dilatation motions of the plate supporting the fuel rods) (Bell and 
Glasstone, 1970). The power oscillations observed in the very compact core of the Canadian nuclear 
reactors (Maple) for the production of radioisotopes, when the steady state power was high enough so 
that the turbulent coolant flow excited strong vibrations in certain structural elements, also seems to 
have a direct relation with the subject of this section. 


Marguet assumes that the period of mechanical oscillations is negligible compared with the mean 
neutron lifetime including delayed neutrons (of the order of 10 s) so that the hypothesis of a constant 
production of delayed neutrons (CDL) mentioned in subsection 2.2, with a suitable version of equation 
(13) may be a good place to start research: 


=f —(1-p)M[9,9.»]-L[4,».¥]+6-M [hy r0.00] oy 


In (67) the field of mechanical variables can be reduced to the scalar volumetric dilation field in the 


solids of the core. If s (7 ) represents the vector displacement field in the solids, then taking its 


divergence we obtain a feedback mechanical variable y is r ) =Ves (¢, x ) to relate the reactivity of 


the reactor with mechanical effects. The fields of thermal-hydraulic variables can be reduced to an 


incremental (relative to the steady state) temperature field w(t, i ) =AT (t, if ) . Through thermoelastic 


effects this temperature field modifies the displacement field in the solids, and besides modifies the 
reactivity acting as a feedback variable. 

A recent simplified mathematical model of a solid spherical fast burst reactor was constructed along 
these guidelines (Kadioglu, Knoll and de Oliveira, 2009). A classic mathematical model due to Bethe 
and Tait, constructed to quantitatively describe accidents in fast neutron reactors, also introduces a 
displacement field as a feedback variable. (Hetrick, 1993) 


In thermal reactors, the diagnosing of mechanical vibrations in reactor core, through the measurement 
of neutron noise, suggests that CDL is a good assumption to study fast nuclear-thermal-mechanical 
vibrations. (Bernitt, 2008)) 

Marguet assumes CDL in the framework of a point kinetic model proposed by Thompson. (Thompson 
and Thompson, 1988). In this section, we first introduce a modification to the model that appears in 
Marguet’s book and then study the possibility, by methods of averaging, of the appearance of nuclear- 
thermo-mechanical oscillations in the framework of the modified model. 

Let us begin with a reactor generating a steady thermal power P, and with the remaining state variables 


at their steady sate values also. Then at the time origin a sudden perturbation is produced in some state 
variable, the reactor power is forced away from its steady state, and a transient begins. 
The introduction of the modified model follows a previous paper (Suarez-Antola, 2007(a)). 
According to the CDL hypothesis, in the equation for power point-kinetics with feedback the term due 
to delayed neutrons appear as a steady source £ - Fy, being B the fraction of delayed neutrons and A 
the mean time between neutron generations (Hetrick, 1993): 

dP 1 


aa ~(p- A) P+©-P, (68) 


The reactivity o will be a function of the temperatures and densities through the corresponding feedback 
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coefficients. However, some temperature effects on reactivity are very fast in comparison with our time 


é : : A ‘ : ; 
scale of reference (this last scale is of the same numerical order as) like certain prompt effects in 


uranium oxide fuels due to an almost instantaneous Doppler broadening in the central region of the 
relatively poor thermal conducting pellets of uranium oxide (Stacey, 2018). The lag between thermal 
power and the temperatures related with prompt effects may be neglected, so that these temperatures 
can be considered as relaxed to equilibrium with the instantaneous power. Temperatures like average 
fuel temperature evolve near the reference time scale, and the lag between them and the instantaneous 
power cannot be neglected. Others, like average coolant temperatures, evolve in scales of much higher 
order, so that these temperatures may be considered as frozen and the heat removal by the coolant may 


be taken as a constant rate /% during the short power excursions studied here (Lewins, 1978). 


The coupling between variations in densities and variations in temperatures necessarily must include 
inertial effects if the time scales of mechanical vibrations in the core structures, t,, , is at least of the 


same order of magnitude as both the time scales £9 of power variations and the characteristic thermal 


time scales t; of those same structures (Boley and Weiner, 1962). If 8 ~ 0.007 and A 10~“ s, then the 
time scale of thermal power variations is = ~10 seconds for a thermal neutron reactor. Vibration 


modes of structural solids like slender bars and thin plates may have frequencies low enough to comply 
with the requirement mentioned above in thermal neutron reactors. The time scale for a fast neutron 


reactor is a thousand times smaller, because in this case A ~10~'s, so higher frequencies structural 
vibration modes could couple also with thermal-nuclear kinetics. 

In the point kinetic model we describe the reactor using the power P(t) , a representative average 
temperature T(t) and a representative average density d(t) of the core, being Py, 7) and dp the steady 
state values of this state variables. 

d(t)— do 


0 


If y(t) = , we suppose that the reactivity is given by: 


p=6p,+a,(P-P,)+a,-(T-T)+a,.y (69) 
The reactivity added from outside, by control bars movement or any other external mechanism, is 69, . 
The feedback reactivity coefficients are: OC; (prompt power) @; (thermal) @, (density). To have static 
stability we assume that both @, and @; are negative and that @ y IS positive. 


For a constant heat removal rate P,, if C =c , -M_ isa suitable thermal capacity (Marguet assumes 
that Mit is the mass of the vibrating component), the representative temperature is given by: 


aT 
C-—=P-B (70) 
dt 
The coupling between dimensionless perturbation of the density of material inside the core and the 
temperature, including inertial effects, is given by the equation of a damped and forced harmonic 
oscillator for the perturbation in the (relative) density: 
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d* d 
Se ee a ee ee (71) 
dt dt 


Here 5 is a thermal expansion coefficient, while c,, is a mechanical damping parameter and @,, is the 


natural frequency of the component which vibrates leading to a change in reactivity. The parameter b 
is a global coefficient of thermal expansion of the component of mass / . All the parameters will be 
considered as constant. 


Besides the excitation of mechanical displacements by thermoelastic effects, the turbulence of the 
coolant flow excites mechanical vibrations in the structural solids of the core (Borsoi, 2001). 


To include this other source of mechanical perturbations a second term could be added to the right-hand 
member of equation (71). This will not be done in the present analysis. 


d eae 
Now let us introduce the new variables: D = a x= if (the so-called logarithmic power). 
t 0 


Define the new functions: 


f (x)=cye* +e," (72a) 
g(x)=e*—-1 (72b) 
The nuclear damping parameters: 
B 
Cy = re (73a) 
(eine a 
er . (73b) 
The thermal-nuclear frequency: 
2 |e; | Lf, 
0, =—— 74a 
"AC a 
The coupling frequency: 
a, .b.P, 
o, =~ (74b) 
: A.C 
From equations (70) and (71) it follows: 
d*v dv o,, .b-P. 
Gn at oe = C *. g(x) (75) 


From equations (69) to (75), eliminating all the state variables with the exception of power, we obtain 
an integral-differential equation (Suarez-Antola, 2007). 
For the purposes of the present work, it is possible to write it as follows: 


d°x dx ; 


ge te) oe eee IK) eC Gn) ane 
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Equation (76) is equivalent to (16) and (17) of subsection 2.4. The impulse response function K (¢ ) is 


of purely mechanical origin and is given by the following equalities, with c, =2-¢,,-@,,: 


K(t)= @'m-G(t) (77a) 
t.. 1 (at _ pet 
G(t)= a (e+ -¢**) (776) 
1 g= Oy {- CaS. -1) (77c) 
[x dt=1 (77d) 


0 


When @, = 0, then O:= 0 In this particular case equation (76) describes a nonlinear oscillator with 


a nonlinear positive damping f (x) and a nonlinear restoring term g(x) , both terms without time 


lags (Suarez-Antola, 2007(a)). 
In this case the steady state is globally asymptotically stable and sustained oscillations of thermal power 
are impossible. 


Ifa, # Oan additional restoring term appears, with a continuous distribution of mechanical time lags 
given by the memory function K (z). Due to the influence of the memory function, the possibility of 
sustained power oscillations cannot be excluded on a priori grounds. A linear analysis in a 


neighborhood of x=0 (P= FP.) with steady power J, as bifurcation parameter shows that if 


CC. 
Bj) ™ 
a Db % b.a, . . 
—— >|a,|.c,, there is a power threshold PR, =~ such that if / > F;, the stationary 
: ee 
b.a 


oa 
operating point of the reactor becomes instable. 


a,,.b 
But if Sar < la, hea then the threshold to instability disappears. 


As shown elsewhere, when the reactor becomes unstable, the power oscillates with increasing amplitude 
until the linear approximation is no longer applicable. (Suarez-Antola, 2007(a)) 

Let us study what happens in this this case applying the method of averaging of subsection 2.4 to the 
nonlinear integral-differential equation (76) of the logarithmic power. 

We begin with the tentative solution for the logarithmic power: 


x(t) ~c(a(t))+a(t)-cosy(t) (78a) 
y(t)=@-t+0(t) (78b) 


The function that gives the offset as a function of the amplitude follows from equation (22). In this case, 
from the definition (72b) of the nonlinear restoring term we have: 


[e(c(a)+a-cosy)-dy = ["(exp[e(a)+a-cosy ]-1)-dy = 0 (79) 
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If I,(z) represents the modified Bessel function of the first kind and order n, the following 


development in series is verified (Abramowitz and Stegun, 1964): 
exp[a-cosy]—1=I,(a)+2: DI, (« -cos(ny) 
Then: 


exp[c(a)+a-cosy |-1=(exp[c(a) |-/,(a)-1)+2-exp[ c( a)|: AC -cos(ny) — (80) 


n=l 


Substituting (80) in (79) we find the function that gives the offset in (78a): 
c(a) =—log,I, (a) (81) 
So: 
g(c(a)+a-cosy)=2 -exp| c(a )]: pyAC -cos(ny) (82) 


n=l 


And: 
f (c(a)+a-cosy) = (cy.e° +c,.e")-I, (a)+ 


+20 7 (83) 
2-¥((-1) Cy.e* to,.e" I, (a)-cos(ny ) 
n=l 
From equations (20a) and (20b) of subsection 2.4: 
[xe )+a-cos(y—a,:t'))-dt'= 
(84) 
2-exp| c(a) |: >, (a )-{K, (n-@,)cos(ny’) + K, (n-@,)sin (ny) 
Here: 
ree o (oe —e 
K.(n-@)= [ K(0)-cos(n-@,t)-dt = - ( . i) (85a) 
0 (w;, - 0) +c -@ 
ies 2 
K,(n-@,)= | K(t)-sin(n-@, -)-dt =» @n (85b) 
0 (@;, -@;) +0 -@ 
The frequency @, is the limit @, = lim,_,, a“ and is determined from the implicit equation: 
@, = @, + @,-K,(@) (86) 
Besides we have the formulae (Abramowitz and Stegun, 1964): 
2 
I,(a)- 1,(a)=— (a) (87a) 
ay i a a’ 
I ie Se ae te 
(a) (5) [ee | ca 


From equations (82) to (87), and from equations (26a) and (26b) of subsection 2.4, with 


ae ede eer eae J K()-8(x-7)-4e 
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we deduce the following equations for the amplitude and phase of the logarithmic power: 


Bs, -£|(6 + Cy B(a))-@? E62) [-20] (88a) 


dt an I,(a 

d I, (a) 

Sy oe 

a o | aut eP) 


Linearizing (88a) the following stability condition is derived for the stability of the reactor in the 
operating steady state: 


>0 (89) 


<a 2-2 (6, iGo (a), Jes (2) FEL) 6(0)|0' (90) 


H(a)=14+-a" + (91a) 
a’ a’ 
—+—— + 
@(a)= 12384 SAG, (916) 
eee 1,(a) 
4 64 


The function IT a) is a growing function of the amplitude, bounded below by 1. 


The function O(a) is a positive decreasing function of the amplitude, bounded above by 1. 
Q, 
If (ae eC; ) = Q, : Blas) becomes negative, the steady state of the reactor becomes unstable and the 
M% 
amplitude of the logarithmic power grows. 
However, its growth is bounded, because in this case there is an equilibrium amplitude @, which is 


stable and corresponds to a stable limit cycle of power oscillation. 


K,(@ 
The existence of this amplitude when (c, + Cr)— @, Ko) <0 stems from the behavior of the 
Op 


2 
C ao, -K(@ 
positive function {Sem EO) ofa) that grows with the amplitude from zero, 
. DM 
d K,(@ d 
being rr positive, until it attains the value Q, De, +¢;) when ae =0 and then 
t M t 


' ’ bias 
continues its growth with hh a negative. 


Equations (81) and (88) define a normal form obtained by the averaging method. 
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When the amplitude a (t ) is small enough (relative tol) Equation (90) can be approximated by the 


canonical equation of the normal form (topologically equivalent to equation (81) and (88a)) that gives 
a supercritical Hopf bifurcation: 


<a(t)=K-a(t)-7-a"(t) (92a) 
Here: 
1 f 2. Kl) _ 92b 
aac = ro) (92b) 
Ea eS) (92c) 


a 
2 M 
From equation (88b) we have an approximate formula for the instantaneous frequency of power 
oscillations, that decreases while the amplitude increases: 


ofa) aon { i lao, f1-<| (93) 


These oscillations are coupled with the corresponding oscillations in the mechanical subsystem, as can 
be seen in Eq. (75), that can be written using the ansatz of the logarithmic power in this way: 

dv... dv...» or bP. 

ton yp ten Da" g(eta-cosy) (94) 
The amplitude of the logarithmic power varies slowly relative to the time scale corresponding to the 
frequency given by equation (86) and to the transient response time of the unforced mechanical 
subsystem. 
Then, once a mechanical perturbation is introduced in a steady state reactor and after a short initial 
transient during which coupled thermal-nuclear and thermal mechanical oscillators adjust their 
dynamics, the mechanical variable will make forced oscillations according to (94). 
Due to the difference in time scale, the speed of offset and amplitude of the logarithmic power can be 
considered as constant while solving (94) for the vibrations of the structural solids. 
In that case the mechanical variable will be dragged by the oscillating reactor power, either towards a 
rest state or towards a stable limit cycle. 


6. Conclusions 


6.1-In the example of a runaway ina simplified model of a nuclear reactor, we saw how to find a formula 
for a threshold amplitude (when it exists) that can be used to determine the minimum amplitude of a 
perturbation that leads to instability, even when there is stability for suitably bounded perturbations. 

A normal form for a static global bifurcation (Equation (38) of section 3) was derived. As far as the 
author knows, it is the first time that it is obtained and applied in nuclear reactor dynamics. 


This example was posed to show as directly and as simply as possible how threshold amplitudes for a 
static instability can be derived from the methods of nonlinear modal analysis developed by Wiktor 
Eckhaus. Because of that it is not intended to be a realistic model to be applied to some practical 
problem. However, the example of the idealized sub-critical reactor with a distributed external neutron 
source, developed in this paper, could be modified and extended to study the space time dynamics of 
sub-critical multiplying systems driven by external neutron sources. This could be a complement, from 
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the restricted nonlinear modal analysis standpoint, to the physics of these systems (Gandini and 
Salvatore, 2001). 


The results of digital simulation reviewed in this paper confirm the dominance of the first mode, perhaps 
unexpectedly far from the bifurcation point. This behavior is related with the existence of an attracting 
slow manifold. 


6.2-In the example related with xenon oscillations we uncoupled the fundamental mode from the first 
harmonic neglecting the coupling terms. 

We reduced the dimension of the state space from six to three. This allowed us to consider in-phase 
xenon oscillations in a mathematical framework quite like the one obtained applying the equations of 
point kinetics with feedback. The analysis done here can be related with the research, using averaging 
methods in a point kinetic model, done by Yoshiro Asahi and Ziya Akcasu. (Asahi and Akcasu, 1973) 
However, in our case, we have a well-defined connection between the static reactivity and the geometry 
of the homogenized core. 


The nonlinear equation for in-phase Xenon oscillations derived by Chernick et al. (1961) is based on a 
model of point kinetics. The nonlinear equation for in-phase Xenon oscillations obtained in this article 
is based on a modal analysis of a spatiotemporal model of the reactor core, so parameters that appear in 
this second case differ from parameters that appear in the first case. 

Equations (57b) and (59a) and (59b) can be considered as normal form equations to study stability and 
bifurcations in the framework of the present mathematical model of in-phase Xenon oscillations. The 
offset term makes them different from the common normal form for a subcritical Hopf bifurcation. 
From (60) we obtained an analytical formula for a threshold of instability to finite perturbations that 
relates the amplitude of a threshold perturbation with the parameters of the homogenized reactor. 


In the second part we uncoupled, in an indirect way, the first harmonic from the fundamental mode. 
Instead of neglecting the coupling terms, we substituted the amplitudes corresponding to the 
fundamental mode by their steady-state values. A three-dimensional state space resulted from this. A 
further reduction in the dimension of the state space was obtained by the assumption that the flux 
amplitude and xenon concentration amplitude are in equilibrium. This gave us a two-dimensional linear 
system equivalent to a damped harmonic oscillator. The damping coefficient and the natural frequency 
was given by a well-defined analytical relationship in terms of the parameters of the reactor. These 
results are a generalization in the framework of modal analysis, of the results for xenon spatial 
oscillations obtained using two-nodes nodal models like the one proposed by Henry in his book. (Henry, 
1975) 

The connection between the assumption about the quasi-equilibrium between the mode amplitudes 
corresponding to the neutron flux and the mode amplitudes corresponding to the Xenon concentration, 
and the theory of slow manifolds, suggested in 4.1 and 4.2, deserves a detailed investigation. 


6.3-In the last example, related with mechanical kinetic effects, the modelling in space-time of the 
previous examples was abandoned by a phenomenological point kinetic model with feedback, after 
some brief comments related with the possible effects of the displacement fields of the solids in the core 
on the reactivity in a homogenized reactor. 

It is reasonable to expect that only the divergence of the displacement field appears in the reactivity. 
Several damped vibration modes will appear in the core, under suitable excitation by thermal-elastic 
effects in the solids or by pressure and shear stress fluctuation effects on the liquid-solid interfaces 
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related with coolant turbulence. A dominant mode of damped vibration could be established but its 
connection with the field of volumetric dilation of the solids in the core is not straightforward. 

So, a lumped parameters model, like the one proposed by Serge Marguet (Marguet, 2917), was 
considered. 


A reactor point-kinetic equation with frozen delayed neutron effects and a mechanical oscillator 
equation for a dimensionless density perturbation driven by thermoelastic effects was constructed and 
studied by averaging methods. 

The original model was modified introducing prompt effects in the reactivity through an additional 
reactivity coefficient that is related directly with the power level in the reactor, without any time lag. 
The nonlinear damping of the resulting nuclear-thermal oscillator is always positive (stabilizing). 
However, when the thermal steady power is greater than a threshold power, the coupling between the 
nuclear oscillator and the mechanical oscillator produces fast (relative to delayed neutrons time scale) 
and growing mechanical and thermal oscillations around an unstable steady state, that tend to a stable 


limit cycle. The threshold of power doesn’t exist if the prompt power feedback coefficient @; verifies 


a,,.b 
the inequality —— <a.) According to the mathematical model, in this case the reactor remains stable 
q y C i g 


m 


at any steady power. 


Equations (88) jointly with the offset formula (81) can be considered as a normal form to study stability 
and bifurcations of the amplitude and phase of the logarithmic power of the reactor. 


When the amplitude a (t ) is small enough, the canonical equation (92) of supercritical Hopf bifurcation 


is obtained. In presence of noisy data, the estimation of parameters in the deterministic model could be 
done following the procedures mentioned at the end of 1.1.2 and 1.2.1. 


A possible generalization of the point kinetic model of nuclear-thermal-mechanical coupling would be 
an elimination of the assumption of a constant production of delayed neutrons (CDL). If CDL is 
eliminated, equation (68) must be substituted by the following: 


BP). ps BF Dla) Pu) 


Here D(t) is the kernel of delayed neutrons (Bell and Glasstone, 1970; Duderstadt and Hamilton, 1976). 


Beginning with this last equation, a generalization of equation of the integral-differential equation (76) 


is obtained. It can be analyzed by the averaging method described in this paper or by methods of digital 
simulation. An application of the averaging method to an analogous integral-differential equation, 
derived to describe BWR power oscillations, can be found in Suarez-Antola (2012). 


Point kinetic models, like the present one, are too crude idealizations for studies in reactor dynamics 
and control, so that detailed numerical and experimental research using full scale system codes for 
specific reactor types must be done to discover the many subtleties that seem to be hidden under the 
umbrella of nuclear-mechanical coupling. Nevertheless, these results could be of some interest to 
discuss aspects of the physical safety of nuclear reactors, related with mechanical vibrations, indeed 
outside the scope of most available neutronic-thermal-hydraulic numerical codes. 
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A necessary condition to use point kinetics in the analysis of a possible thermo-mechanical coupling 
due to thermal-elastic effects in light water reactors is a homogenization of a highly heterogeneous 
reactor core. From a neutron physics point of view, this is not unreasonable due to the long mean free 
paths of neutrons in the reactor core relative to the diameters of and the distance between fuel pins in 
fuel assemblies. The main problem with the present model is related with how the thermo-elastic effect 
and its consequences appear in the core. The process begins with a temperature change in the fuel 
pellets. If this change happens in a time scale of an order of magnitude (let us say, 107 s or less) less 
than the order of the time scale of the heat transport from fuel to coolant (let us say, 10°! s), a coupling 
between the temperatures field and the densities field in core materials like the one assumed in the 
present model seems highly improbable. 

However, if the average temperature in the fuel changes in a time scale of 10°! s, and this change appears 
jointly with low frequency mechanical vibrations that introduce long range mechanical correlations, 
and that can be excited by distributed thermo-elastic effects in several fuel elements, then a point kinetic 
model, as a first approximation, could be acceptable if the mechanical oscillator is related with a 
fundamental mode of vibration. In any case, an in-depth theoretical study of the possibility of nuclear- 
thermal-mechanical coupling must be done using suitable systems of nonlinear partial differential 
equations, adapted for each type of reactor. 
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